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On the Application of Statistical Concepts 
to the Buffeting Problem’ 


H. W. LIEPMANN?T 
California Institute of Technology 


(1) ABSTRACT 


This study is intended as a step in the direction of an under 
standing of the problem of buffeting. The buffeting discussed 
here is of the type encountered when the wake of some part of 
the airplane interferes with the tail surface, as happens, for 
example, during the stall. To approach the general buffeting 
problem immediately does not seem promising; hence, a simpler 
problem is studied to bring out some general features of the 
buffeting problem which are believed to be important. The 
problem studied is the lift force exerted upon a two-dimensional 
thin airfoil moving in turbulent air. The state of turbulence 
will later on be restricted to isotropic turbulence, since the 
characteristics of this type are the best known and simplest. 
The explicit results gained from the present analysis will not be 
immediately applicable to tail buffeting, since the model is cer- 
tainly oversimplified. 

The lift exerted upon an airfoil moving in turbulent air is re 
garded here as the response of a linear system to a random 
forcing function. The explicit use of the random character of 
the forces exerted by the turbulent air is believed essential for 
an understanding and handling of buffeting. This concept 
obviously extends beyond the applicability of the simple model 
discussed here in detail. The concepts used in studying the 
Tesponse of a linear system under the action of a random input 
force are by no means new. For example, the idea of random 
phases which leads up to the so-called power spectrum goes back 
toLord Rayleigh.! The theory of generalized harmonic analysis 
has been developed in connection with communication engi 
neering and related problems, mostly by Wiener® * and his school. 
A parallel and often independent development went on in the 
theory of turbulence, where the same concepts of correlation 
functions, power spectrum, etc., have been systematically intro 
duced by G. I. Taylor,4 von Karman, and others and used 
&xtensively since then. In other branches of aeronautics and 


Presented at the Aeroelasticity Session, Twentieth Annual 
Meeting, I.A.S., New York, January 28-February 1, 1952. 

*H. Luskin introduced the author to the problem of buffeting. 
The approach presented here owes much to discussions with him. 
Thanks are also due to Drs. J. D. Cole, P. A. Lagerstrom, and 
C. DePrima for constructive criticism 

} Professor of Aeronautics, Guggenheim Aeronautical Labora- 
tory, California Institute of Technology. Also Consultant, 
Douglas Aircraft Company, Inc., Santa Monica. 


aeroelasticity, these concepts and their use are comparatively 
unknown, and this fact excuses the rather lengthy and detailed 
discussion of simple concepts in the present paper.t The author 
is well aware of the fact that many of the problems discussed 
and many of the concepts used are not new. It is the main 
purpose of the paper to call attention to the use of the statistical 
concepts in aeronautics and to demonstrate their application 
with a simple example. For a further study of these and similar 
subjects, as well as for a more complete list of references, the 


reader is referred to Wiener’s work* and a paper by S. O. Rice.® 


(Il) THE RESPONSE OF A MECHANICAL SYSTEM 
TO A RANDOM FORCE 


aoe DEVELOP THE GENERAL IDEA underlying the 
buffeting analysis, it is most convenient to look 
first at a simple problem that has been treated before 

namely, the motion of a pendulum in a turbulent 
The same type of analysis covers many more 
for example, the response of an 


stream.’ 
similar problems 
acoustical resonator to random noise, the response of 
an electrical circuit to electrical noise, Brownian motion 
of a torsional balance, etc. These problems all have in 
common the following: The given system is described 
by a simple differential equation of the oscillation type. 
That is, let y(t) represent the angle of the pendulum as 
function of time, the displacement of an oscillating 
mass, or the current in an electrical circuit. y(t) is 
determined by the differential equation 


d*y dy 
+ 


f F(t) (1) 


dt? d 


where £8 represents a damping parameter, wy 
eigen frequency of the system, and F(t) is an impressed 


+ wy -_ 
t 


is the 


t Reference should be made to two unpublished M.1I.T. 
theses!* 13 in which the problem of gust response of an airplane 
and of atmospheric turbulence was attacked using statistical 
concepts similar to the ones used in the present study. The 
author became aware of these theses only recently. For the 


engineering approach to buffeting see reference 14. 
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force. If F(t) is given explicitly as a function of time, valid for the respective power spectra 
we deal with a simple inhomogeneous differential equa- F(a) f( 

Ww w) 


In the simplest case let 


F(t) = Re{ A-e'*'} 


tion. 


that is, the force is a periodic function of ¢. Then the 


steady-state solution y(t) is found by writing 
y(t) = Ref B-e'“'t! (A, B real) 
B and ¢ are then determined from Eq. (1) 
tan g = — Bw/(an”? —w?) (2) 


, A* ' 
aiiien hr pene? reer ee (3) 
(wo — w*)* + B*w* 
The mean square of y(t) is related to the mean square 
of F(t) by Eq. (3), since 


y? = B?/2; fF? = A?/2 
and, hence, 


7 fF? F? 
Pe ae ae ; (4) 
(wo” — w*)* + B*w* |2(w)|* 
lo - e e . 

|z(w)|? is the absolute square of the impedance of the 
system. 

Now assume that F(t) is not a well-behaved sine 

wave and not periodic, but a random function. That 

is, F(t) is defined only by certain mean values— for 


example, F?, F, etc., where 


” 1 T 
f* = oT fs. F(t)” dt 


with 7° sufficiently large. The actual trace F(t) will 
vary in an irregular fashion, and, hence, it becomes 
meaningless to ask for the trace y(t). Especially, the 
phase of y(t) will lose all its meaning. 
question, however, remains: How is the mean square 
of y related to the mean square of F? In the case of 
the pendulum in turbulent air the question is: What 
is the mean square deflection of the pendulum if the 
mean square velocity fluctuation (the turbulence level) 
is known? 

From Eq. (4) it is evident that, in order to answer 
this question, we have to know some of the frequency 
properties of F(t). The necessary information is fur- 
nished by the power spectrum f(w). f(w) is defined in 
the Appendix; f(w) satisfies the relations 


F2 = f fo)do « P= f f(w) dw (5) 
0 0 


that is, the first of Eqs. (5) represents the contribution 


The following 


to F* arising from frequencies less than w. Similarly, 


we may define a power spectrum of y by 


y= f g(w) dw (6) 
0 


Generalized Fourier analysis (for example, references 2 
and 3) states that, for the case of a function defined by 
its power spectrum and mean values, Eq. (4) remains 


(oe + - z 
(wo — w*)* + B*w° Z(w)| 


Hence, from Eq. (6) 


: f(w) dw 
y~ = : (X 
0 |3(w)|° 


Eq. (8) includes Eq. (4) as a special case. 

The relation, Eq. (7), between the power spectra oj 
the output g(w) and the input f(w) is true for any linear 
system. 2(w) denotes the impedance of the linear 
system, which, of course, is not necessarily of such a 
simple form as Eq. (4) or (7). 

Consequently, in order to obtain the response of a 
system to the random forcing function, one has to 
know the impedance of the system and the power 
spectrum of the forcing function. Then the mean 
square of the response is obtained from Eq. (8). A 
proof of this relation is included in the Appendix under 
heading (6). 

If f(w) is a slowly varying function compared to the 
variation of the impedance with w, Eqs. (7) and (s 
can be written approximately as 


—— 7 ; dw . 1 f(a) 
y? = f(wo) - a ; «9 
0 (w- — ei)" + B-w- 4 Bar- 


Hence, the mean square response depends upon the 
value of the power spectrum of the force at wo, the 
eigen-frequency of the oscillator, and the damping 
characteristic of the system. Eq. (9) can be thought 
of as defining f(w) by a measurement with a wave 
analyzer. It is a known and useful equation, for 
example, in the theory of Brownian motion. 


(III) MorTion oF A Two-DIMENSIONAL AIRFOIL 
IN INCOMPRESSIBLE TURBULENT FLOW 


Consider a flat thin airfoil of chord c moving with a 
constant velocity U through turbulent air. Let x lie 
along the chord, y along the span of the wing, and z 
normal to span and chord. The components of the 
fluctuating turbulent velocity u, v, w are assumed small 
compared to U. Because of these turbulent fluctua- 
tions, a time-dependent apparent angle of attack a 
exists at the airfoil, and, hence, fluctuating lift forces 
are produced. The fluctuating angle of attack a is 
given by 


a=w/U (10 


as long as the fluctuations are small. a(t) now plays 
the role of the forcing function. The “‘response’’ 1s 
the fluctuating lift of the airfoil or, better, the lift 
coefficient C,(t) 

In an approximate way the considerations of Section 
(II) can be carried over to the present problem. To do 
this it is first necessary to define an “impedance” of 


the airfoil. Sears® has studied the motion of an airfoil 
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APPLICATION OF 
to a sinusoidal gust, and his results can be used 
for our purpose. Sears considers the fluctuation of the 
velocity of the form 
w(x jt)/U — an ett -@ U)] 

that is, the vertical velocity consists of harmonic waves 
of angular frequency w and wave length 27U/w—.e., 
a fixed pattern of sinusoidal gusts is passing the airfoil 
at a velocity U, the flight velocity of the wing. For 
this velocity distribution Sears finds the following ex- 
pression for the lift coefficient: 


Cz) = Qray-e'y(k) (11) 


where k = we/2U is essentially the ratio of chord 
length to wave length of the gust. Thus, the function 
y(k) has the nature of an admittance, relating the forc- 
ing function a(t) to the response C(t). 

Turbulent fluctuations are essentially three-dimen- 
sional—that is, u, v, w will be functions of x, y, 2, ¢. 
For a complete analysis one would therefore require a 
three-dimensional formula equivalent to Eq. (8) and a 
result of three-dimensional airfoil theory equivalent to 
Eq. (11). For the present it seems sufficient to con- 
sider only the component w and the dependence upon 
x,t. Thus, in turbulent flow we consider a fluctuating 
velocity or angle of attack of the form 


a(x, t) = w(x, t)/U 


If it is now assumed that the turbulence pattern does 
not change appreciably during a time of the order c/ U, 
then the turbulent angle of attack will also depend 
upon ¢ — (x/U) only and Sears’ results can be applied. 
This assumption is frequently made in turbulence 
analysis and requires essentially the condition that the 
rate of change of the velocity following the particle is 
small compared with the rate of change of the velocity 
at a fixed location. 

With this assumption we can carry over the analysis 
of Section (II) directly. If the power spectrum of a(t) 
is denoted by f(w), we can thus write, for the mean 
square lift coefficient, 


Cr? = 4r’ f f(w)| e(k)|? dw 
0 


The function g(k) is rather complicated 
2 Jo(k) Ki (tk) + iJ (k) Kot(k) 
7 Ki (tk) + Ko(tk) 


¢g(k) (12) 
J, and K, are Bessel functions. 

For the present analysis we require \y|?. To obtain 
¢|* from Eq. (12) is fairly involved. The resulting 
expression is complicated, and it seems more appro 
priate for the present study to find a simple approxi- 
mation to |g|2. For this purpose |¢\? was plotted from 
Sears’ paper as a function of k. Fig. 1 demonstrates 
that |y|? is well approximated by 


||? = 1/(1 + 2xk) (13) 
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x 
Fic. 1. 


Eq. (13) results as an interpolation formula: | ¢(0)|? = 1 
and using the asymptotic expression in Eq. (12) 
\g(k)|* > 1/2rk for large k. Using this function the 


mean square lift coefficient can be written 
f(w) dw 


+ (#e 


(14) 


U)w 


f(w) denotes the power spectrum of the cross component 


that is, 


— w? ia ; 
a =—_—= f(w) dw 
U?? 0 : 


Now for isotropic turbulence it has been found (for 
example, references 9 and 10) that the bulk of the 
turbulent energy is contained in a spectrum region that 
can be described by the simple power spectrum 


of turbulence 


, w? L 1+ 3(L%w?/U?) 
fiw) =->- rm an 


: wert (15) 
aU {1 + (L2w?/U?)]? 


where L denotes the so-called ‘“‘scale of turbulence’”’ 
[see the Appendix under heading (2) for definition]. 
Call Lw/U = &; then, 
f(w) = f(0) (C1 + 3é)/(1 + &)?] 
Using Eq. (15), Eq. (14) can be integrated explicitly. 
The result is 
w? in — fF n+3 | 
C,? = 44° — + (nlogn? + x 
; a Peay 2r(n? + 1) — 
(16) 


with rc/L = 1. 
coefficient of a two-dimensional airfoil of chord ¢ in 


Eq. (16) gives the mean square lift 


isotropic turbulence characterized by the intensity 


w?/U? and the scale L* (Fig. 2). 

* Actually, one also has to take into account the fact that the 
turbulence is not two-dimensional. That is, there also exists an 
effect due to the fact that different parts of the airfoil along the 
span will find themselves in different phases, and, hence, the 
overall result is less than Eq. (16). The effect is similar to the 
length correction of a hot-wire anemometer. 
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Fic. 2. 


Evidently, if c/L ~0, we have an airfoil of small 
chord in large-scale turbulence and 


° 
7p)- 


ew 
C12 > 49? — = 42a? 


The airfoil behaves quasi-stationary. If, on the other 
hand, c/L becomes large, we deal with an airfoil of long 


chord in small-scale turbulence. It follows from Eq. 


(16) that C,*—> 0. 
each other completely and the net lift is zero as might 
be expected. 

A phenomenon not included in an analysis using a 
spectrum of the form (15) but possibly of major impor- 


That is to say the “gusts” cancel 


tance for buffeting is the so-called ‘intermittency’ in 
wake flow.* That is, the edge of any wake fluctuates 
with a large-scale motion so that a point situated near 
the edge will sometimes be within the wake and some- 
times outside the wake. For a tail surface situated near 
the edge of the wake of a stalled or partially stalled 
wing, this “‘intermittency’’ may thus be extremely im- 
portant in determining the lift forces. A detailed 
analysis would require much more experimental study 
than is available at present. For a crude idea of the 
effect, consider the flow at the tail as a region of uni- 
form downwash switched on and off at irregular inter- 


vals. Such a flow is probably a good model for the 
conditions in the wake of an intermittently stalling 
wing. If the probability of switching over from one 


régime to the other is assumed to be governed by a 
Poisson distribution, one has again a simple form of the 
power spectrum 


w? 2L I 


U2 eU 1+ (wL/U)? aa 


f(w) 


where w denotes the downwash velocity and L/U is 


the average length of time the downwash is ‘‘switched 


on.’’ The formula corresponding to Eq. (16) is then 


P w? 2 n log n + (2/2) 


C.? = 4x? — 16a) 
L U2 x ve: (l6a 





* For the concept of “intermittency” see, for example, Town- 
send, A. A., The Eddy Viscosity in Turbulent Shear Flow, Phil. 
Mag., Vol. 41, p. 320, 1950. 


SCIENCES—DECEMBER, 


(IV) AEROELASTIC EFFECTS 


If elastic effects are to be considered, one may ask 
for the deflection of the wing under the lift forces 
Again an ‘‘impedance’’ due to the elastic structure of 
the wing will enter. Assume, for example, that the 
elastic character is of the simple type as expressed jn 
Eq. (1). The impedance corresponding to this type of 
structure has the form 


(w? — a”)? + Bw? (17 


and, hence, the mean square deflection y*, due to the 
fluctuating lift, will be obtained from 


A = [ : f(w) dw 
yo =4r ; i: ae ——= oe 
0 [1 + r(c/U)w][w? — wo)? + B2w?] 


In more general terms we have 


ym f f(w)¥(w)x(w) dw (19 
0 


where f(w) denotes the power spectrum of the turbu- 
lence, ¥(w) is the absolute square of the admittance of 
the wing, and x(w) is absolute square of the admittance 
of the elastic structure. 

deflections will, of introduce 


Structural course, 


additional aerodynamic effects. For example, oscilla- 
tion of the fuselage due to the tail loads will introduce 
plunging motions of the tail surface and thus aero- 
dynamic damping of the oscillation. These coupling 
and feedback effects will vary from case to case and 
may sometimes be rather complex. 

From Eq. (18), for small damping, ‘follows the simple 
relation: 


? . f(a) 
yr = 4r° —— : (20 
2Bor? [1 + a(anc/ LU) ] 


For the simple case of the input spectruin (15) this yields 
9 l + 3(L7 a? U?) (21 


U [1+ a(anc/U)][1 + (L2a?/ U2) |? 


A result like Eq. (20) applies within the limits of the 
analysis, for example, to the response of an elastically 


restrained flat plate used for turbulence measurements.] 


(V) RESPONSE TO A VORTEX STREET 


Sometimes buffeting is assumed to be caused by a 
von Karman vortex street. While this may be the case 
in specific instances, it does not seem to be the general 
case. Even if a vortex street of nearly regular frequency 
is present, little will be gained by a detailed analysis 
involving the phase of the lift, as well as its magnitude 
From the general formulation for C,? given here, the 
mean square lift due to a vortex street is, of course, 


obtainable. The velocity components in a_ regular 


+ A device similar to this has been used by Ackeret for turbu 
lence measurements in water (verbal communication ) 
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APPLICATION OF 


yortex street are sinusoidal, and, hence, they have a 


power spectrum / 

f(w) = (A?/2)6@~ — Q) (6=6 function) (22) 

where 2 represents the vortex frequency. Inserting 
this into Eq. (14) yields 

Ci? = 29A*{1/[1 + (mc/U)Q)} (23) 


This result could have been obtained directly from 
Sears’ analysis. 

In most cases even a vortex street shows a continuous 
background spectrum. This is clearly seen in Fig. 3, 
where the power spectrum as measured behind a cylin- 
der at a Reynolds Number of 4 X 10° is shown. Note 
the distribution of energy between the continuous part 
and the discrete band. The spectrum shown is the one 
corresponding to the fluctuation in the direction of the 
mean flow. Hence, the simple analytical spectrum 
corresponding to Eq. (15) is simply of the form 


f(w) = a/(1 + bw?) 


compare with the Appendix under heading (5)]. The 
relative energy in the shedding frequency band and 
the continuous part is proportional to the relative area 
as indicated on Fig. 3. 

The continuous part of the spectrum obviously be 
comes more predominant further downstream from the 


cylinder. * 


(VI) THE PROBABILITY AND FREQUENCY OF 
EXCESSIVE VALUES OF THE RESPONSE 


It will often become necessary to ask not only for the 
mean square response but also for the probability that 
This question is 
If the probability 


a certain value of y is exceeded. 
natural for strength considerations. 
distribution ot the response is known, then this answer 
can be given immediately. Often it will be possible to 
assume, at least approximately, a Gaussian distribution 
of the response. In this case the probability that the 
variable y will exceed k times the root mean square 
value V y?, for large k, is 


PG 2) (k2) 


kV 20 


For example, the probability to exceed 3V y? is about 
0.002. One can also give some general, but also broad, 
limits for the probability that a certain value is exceeded 
without knowing the distribution explicitly. For ex- 
ample, according to Chebycheff’s inequality,'! one has 


prob. y > RV yt ~ (24) 


prob. ‘ y-y> RV y2 I< 1/k? (25) 


* Fig. 3 was taken from a paper “On the Development of a 
Turbulent Wake from a Vortex Street,’’ by A. Roshko, prepared 
under N.A.C.A. Contract NAw-6048 at the Guggenheim Aero- 
nautical Laboratory, California Institute of Technology. This 
paper will appear soon as an N.A.C.A. publication. 
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The probability that y differs from its mean value by k 
times the root mean square is less or equal to 1/k?. 
Applying this to the lift coefficient for the wing as 


computed here yields 
prob. ' Ey > eV ci3t < l k? 


Thus, C, will be ten times larger than VC)? in less 
than one out of a hundred cases. The Chebycheff in- 
equality is known to be far too broad for most prac- 
tical applications, and the upper limit given is, in gen- 
eral, much too high. The estimate can be considerably 
improved if one knows more about the distribution 
function. In the present application, for example, it is 
safe to use Gauss’ estimate with zero skewness, !* which 
gives 


(26) 


prob. ‘ Cr = eV Cu < 4/9R? 


Further improvements on this estimate are possible." 


For fatigue considerations it will also be useful to 
ask for the number of times per unit time a certain 
In a number of cases the answer to 
The average 


value is exceeded. 
this question can be given explicitly. 
number of ‘‘¢ values’’ NV; of a stochastic function /(t) 
that is, the number of times /(/) passes the value & 


+ « 
N; = [ p(é, )\n| dn (27) 


- 


is given by ® 


where p(£, 7) d§dn denotes the probability of finding 
I(t) between ~ and ~ + dé with a derivative /’(t) be- 
tween 7 and 7+ dn. If J(t),/’(t) are Gaussian, N; 


becomes simply 


N; = (1, ne® 262 oe 
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The average number of times the function exceeds 


nV 2 is thus given by 


‘NN = Love [a (28) 


The factor '/2 stems, evidently, from the fact that J(t) 
passes the value é half of the time with positive slope 
and half of the time with negative slope. Only the 
In terms of the power spectrum 


y| g(a) do] 
0 


(29) 


former enters here. 
g(w) of J, the result is 


— k2 9 
1 


" Make peer 
N = | w'g(w) dw 
0 


» 


“ kv/ 2 27 


(VIL) CONCLUDING REMARKS 


It is believed that, in spite of the simplicity of the 
model discussed here, a few important features of the 
buffeting problem have been brought out as follows: 

(1) Since one deals with a phenomenon that is 
largely random in character, there is no use in making 
an analysis that will give the phase of the motion. 
Only mean quantities are defined and useful. 

(2) The use of the power spectrum concept for the 
motion of the air and of the impedance concepts for 
the aerodynamic and elastic response makes the 
analysis relatively simple. 

For a complete theory of buffeting, obviously much 
more analysis is required and, especially, considerably 
more knowledge of the turbulence structure of the flow 
in wakes. The latter is very much an experimental 
program. However, it is likely that, once there are 
some data obtained on the flow in the wake of a stalled 
wing, for example, one may find simple rules for the 
power spectrum and turbulent intensities. For ex- 
ample, the turbulent intensity must be related to the 
drag of the body that produces the wake; the scale, 
to its lateral dimensions. 

It should be pointed out that there exists a number 
of other problems in aeronautics for which the concepts 
used here are the natural ones. For example, wave 
drag of rough surfaces and, in general, the definition 
and measurement of aerodynamical roughness belong 
here. The motion of smoke stacks and steel towers, as 
well as airplanes in gusty air, also belongs to this type 
of phenomenon. (The gust problem for aircraft has 
recently been attacked using the technique of station- 
ary random processes.’ '*) Taxiing over rough ground 
is another problem in the same class. 


APPENDIX 


(1) Mean Values 

Assume that velocity fluctuations in turbulent flow 
are observed with a hot-wire anemometer and are 
reproduced on the screen of an oscilloscope. Let u(f) 
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be the fluctuation measured at a particular point in 
the flow. u(t), for example, as viewed on the screen of 
the oscilloscope, will show an irregular character, and 
no periodicities will be apparent. If an instantaneoys 
picture of the screen is taken repeatedly in intervals, 
one will find that the trace u(t) will differ from picture 
to picture. Hence, it does not make sense to ask for 
the complete function u(t); only certain mean valye; 
are reproducible and measurable quantities. One may 
find the mean values of u(t) in the following way 


Take a large number, say JN, pictures of the screen and 


note the value of u(t) at the center of each picture, | 


This gives N values for u; The mean value 7 is given 
by 


u(t) = (uw + ue+...uy)/N 


To be more precise one should define 7(t) as the limiting 
. 5 
value for N > o, that is 


N 
a(t) = lim (>ou,/N) 
0 


N— 
For turbulent flow u(t) will be negative as often as 
positive and a(t) = 0. However, the mean value of 
u*(t) will, in general, be finite and represents a measure 
for the intensity of turbulence or at least of the par- 
ticular component of turbulence measured. Again, 


u*(t) is defined by 


u*(t) = lim (dou? 'N) 


N= « 
Now, if turbulence is observed at a fixed position in a 
steady mean airflow, then the fluctuating field will be 
stationary in the sense that the mean values do not 
depend upon the time—that is, uw? at a specific point 
in the field will have the same value at any time 
Clearly, many flows of interest have this property; 
the flow in the wake of a stalled airfoil obviously has 
this property if the stall persists for a large length 
of time compared with the characteristic times of the 
turbulence in the wake. This will certainly be true for 
a slow stall and even in many cases of pull-out. 
For the analysis in this paper it is always assumed 
that the fluctuations are stationary. 
that in this case the mean values defined above in the 


It can be shown 


way of ensemble averages are identical with the average 
taken over time provided the time scale T [see (2) be 
low] is finite, and similar integrability conditions apply 
to higher order correlations. 


(2) Correlation Coefficient 


Assume again that N pictures of the turbulent trace 
u(t) are taken. Take the value of uw at a particular 
point on each picture (for example, the center) and call 
it u(t). Then take the value of uw at a neighboring 
point ¢+ 7, say, and call this w(¢+ 7). Form the 


average of the product of the two—that is, 
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APPLICATION OF 


\ 
u(t) u(t + r) = lim [Dou,(t) u(t + 7)/N] 


N—> « 0 


be independent of ¢ and a function of + only. This 
means that the mean product will depend only upon 
the distance between the points chosen on each picture 


but not where they are chosen. Hence, we find a func- 


tion 
N 
¥(r) = lim [Sou,(t) u(t + 7)/N 
V— 0 


¥(r) is called the correlation function or, more exactly, 


autocorrelation function. 


It is easy to establish some properties of y. If r = 0, 
we have 
V N 
du; (t) u(t + 7) your (t) 
¥(0) = lim ” : =lim ~— =u" 
\— N N— « N 


If r becomes large, we will often expect y — 0, since, 
if the process under consideration is random in char- 
acter, we do not expect any correlation between values 
of u(t) at tand, say, of w(t) a half an hour later. Hence, 
if r is large, it is equally likely that a positive value of 
u(t) will be multiplied by a positive or negative value 
u(t+ 7), and, hence, the sum will come out to zero. 
Hence, we expect ¥(7) to have a behavior as shown in 
Diagram A. Whether the function approaches zero 
monotonically or oscillating depends upon the specific 
Furthermore y must be an even function 

Clearly, ¥(7) is some measure 
that is, the time distance 


process. 
that is, ¥(r) = ¥(—7). 
of the ‘‘scale’’ of the process 
over which there exists a dependency between mean 
In this sense y has a mean- 
To define the 


values of the fluctuation. 
ing for the definition of ‘“‘eddy size,”’ etc. 
“scale,’’ or better ‘‘time scale,’’ of the motion more 
precisely, one defines a time 7° by the area under the 


correlation curve and writes 


¥(0)T = Fj V(r) dr 


(3) Space and Time Correlations 


So far, the correlation has been defined in time only. 
Clearly, in homogeneous turbulence one may also de- 
fine stationary correlations in space. Instead of choos- 
ing ¢ and t+ 7, one may measure u(t) at x, y, 2 and 
«+ Ax, y + Ay, 2 + Az and form a correlation func- 
tion; or one may pick a different velocity component 
v(t) or w(t), and, hence, several different correlation func- 
It is here that the simplification 


tions may be defined. 
Here, as 


due to an isotropic turbulent field enters. 
shown first by Taylor and by von Karman, two corre- 
lation functions are sufficient to give all of the others. 
The two chosen are ordinarily denoted by f and g and 
are defined in easily understood picture language by 


Diagram B—that is, f(r) is the correlation between 
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DIAGRAM B. 


velocity components along the line in the direction of the 
component. g(r) is the correlation between points 
chosen at a distance r in the direction normal to the 
We have 


direction of the velocity components. 


u(x, y,2,t) ulx+r,y,2,t) = f(r) 
u(x, y,2,t) uUlx,y + 7, 2, t) g(r) 


v(x, y,2,t) vx+r,y,2,t) = g(r) 

If the field of fluctuation u, v, w is superimposed upon 
a mean flow of velocity U in the x-direction and, fur- 
thermore, if u*/U*, v?/ U*, w?/U* are small quantities, 
then it is oftea possible to interchange time and space 
variables and consider the field of turbulence simply 
transported along with the velocity U along x. Conse- 
quently, the time 7 entering into the correlation func- 
tion can be replaced by r/U, where r is chosen along 
If this is done, there results a unique rela- 
For 


the x-axis. 
tion between (7) and g(r) and f(r), respectively. 


example, 
u(t) u(t+ 7) = (rt) = f(r/U) 
v(t) v(t+ 7) = (7) = g(r/U) 


(4) Concept of Power Spectrum 


Assume that the time correlation function y(r) is 
Form the cosine transform f(w) of ¥(7), 


») » 
| ¥(7) cos wr dr 
TJ0 


This function f(w) is identical with the power spectrum 
This becomes plausible if one ob- 
For 


known. 


f(w) = 


introduced before. 
serves a few properties of ¥(7) as defined above. 
example, using the inversion formula for a cosine trans- 


form, we have 


¥(7r) = f f(w) cos wr dw 
0 


and, hence, 
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¥(0) = f f(w) dw 
0 
Since 
¥(0) = uv? 
we have 


u? = i f(w) dw 
0 


(5) Simple Forms of the Power Spectrum 

Experiments® have shown that in isotropic turbu- 
lence the correlation functions f(7) and g(r) have nearly 
the form 


f(O)e~"’” 1 
g(O)e~”’” [1 — (r/2L)]f 


f(r) 
g(r) 


r> 6 


Hence, the corresponding time correlation functions are 
ti(r) = (Oe 78 
o(r) = yo(O)e~74/* [1 — (U/2L)r] 


The corresponding spectra are 


I 
fi(w) al fi(0) | + & : _ wl 
1+ 32 ~ U 
»(w) = fo(0) —— 
folw) = fol (+ 2)? 


In the analysis we dealt with the angle of attack dis- 
tribution in the x-direction [Section (III)]. In this 
case we want 


w(x, t) we +7, 0) 
that is, g(r) and, hence, Yo(r) or fo(w). 


(6) Proof of Eq. (8) 


For completeness, a short proof of Eq. (8) will be 
is convenient here to write the 


Let 


included here. It 
Fourier integrals in the complex form. 


L(y) = F(t) 


denote a linear differential equation relating y(t) to a 
forcing function F(t). Denote by f(w), g(r), the power 
spectrum and autocorrelation function of F(t); and by 
g(w), ¥(7), the power spectrum and correlation function 


of y(t). We have 


Fi = [ f(w) dw = ¢(0) 
0 
y? = f g(w) dw 
0 
iu 
er) =5 f f(w)-e*” dw 
I aie 
¥(7r) = if g(w)-e'*" dw 


¥(0) 
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Call h(t) the response of the linear system to a unit 
impulse at ¢ = 0. 
and for a system with positive damping, the solution 
y(t) can be written 


t 
y(t) = / F(r) h(t — r) dr 


Since h(t) is zero for t < 0, the upper limit of the integra] 
can be extended to +-. 
changed to o, say. Thus, 


+« 
y(t) = J F(r)h ({-—7r)dr = 
+ oO 
f F(t — o) h(c) de 


W(E) = y(t) v(t + &) 


For a process going on since { = ~ q 


The variable ¢ — + can be 


The correlation function 


using the above integral becomes 


y(t) y(t + £) = 


ff e 
or 
4 «x 
y(é) = JS g(— + o — a’) h(c) h(o’) dada’ 


Replacing y and ¢ by the Fourier integral involving the 


F(t — o) F(t + & — a’) h(c) h(o’) da da' 


respective power spectra yields 


f g(w)e'* dé = 
+ « P 
SSS f(w)e'* h(a) -e*” -h(o’)e—**”' da do’ 


Now the impedance Z(w) is defined as the ratio of input 


over output for a sinusoidal input. Hence, 


l _ 
(wn) = f h(a)e’’’ da 
2(w —« 
twt si iwt l l 
g(wje* = f(w)e" dw 
— « 2(w) 2(w)* 


9 


and thus 


J 
g(w) = f(w)/|z 
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Calculation of the Stability of the Laminar 


Boundary Layer in a Compres 


sible Fluid on a 


Flat Plate with Heat Transfer’ 


E. R. VAN DRIESTT 


North American 


ABSTRACT 


The amount of cooling required to stabilize completely the 
laminar boundary layer on a flat plate at supersonic speeds with 
zero pressure gradient is calculated using the Lees-Lin theory. 
It is found that, above about Mach Number 9, it is impossible to 
stabilize the layer with any amount of cooling when a Prandtl 
Number of 0.75 and the Sutherland viscosity-temperature law 
The distribution of mean properties across the 
The inviscid, 


are assumed. 
boundary layer are obtained by the Crocco method. 
as well as viscous, solutions are determined. Use of the stability 
results in missile aerodynamic heating analysis is illustrated by 
means of V-2 flight-temperature data. 
NOMENCLATURE 

Variables, Parameters, and Functions 

x,y = geometrical coordinates parallel and perpendicular 

to the flat plate, respectively 
velocity components of the boundary-layer flow in the 


“v0 = 
x and y directions, respectively 

p = mass density 

m = coefficient of viscosity 

k = coefficient of thermal conductivity 

Cp = specific heat at constant pressure 

1 = enthalpy 

r = absolute temperature 

T = shear stress 

I,J = functions connected with enthalpy distribution 

S = Sutherland constant, 198°R. 

6 = ratio of Sutherland constant to free-stream temper- 
ature = S/T, = 0.505 

Y = ratio of specific heat at constant pressure to specific 
heat at constant volume 

Cr = wave velocity of disturbance 

ac; = amplification factor 

p = pressure 

v = kinematic viscosity 

Pr = Prandtl Number 

M = Mach Number 

Re = Reynolds Number 

Re, = pud*/p 

Reza = Potta%/a 

Re... = (Rerq)' 

q = amplitude of disturbance 

ny = wave length of infinitesimal disturbance 


a = 27/X 


Presented at the Aerodynamics Session, Twentieth Annual 
Meeting, I.A.S., New York, January 28-February 1, 1952. Re- 
vised and received June 27, 1952. 

* The paper is essentially AL-1334, North American Aviation, 
Inc., under the same title. The author gratefully acknowledges 
the assistance of Mrs. Phyllis Watte’, Mrs. Betty Pfouts, and 
Mrs. Doris Evans in carrying out the computations. 

t Research Specialist, Aerophysics Laboratory. 
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Aviation, Ince. 





ay = 2rL/x 

a = 276* 

Se = 2V'x/pottot n't 

Vx = v/Ve 

Uy = u/Ue 

Ce = Cr/Ug 

Py = p/Pa 

T= T/T 

t. = 4/1. 

My = U/Mo 

L = x/V potlaX/pe 

Ve = normal distance parameter, (y x)V px Up X/ Me 

2 

5* = displacement thickness f, (1 — p,u,) dy 
Subscripts 

00 = free-stream conditions 

W = wall conditions 

1 = based on displacement thickness 


INTRODUCTION AND SUMMARY 


4 nw QUESTION OF WHETHER the boundary layer of 
an aircraft is laminar or turbulent is all the more 
important today because of the high friction tempera- 
tures encountered at supersonic speeds and consequent 
effect on the structural design of the vehicle. Since 
turbulent flow results in considerably greater skin- 
friction and heat-transfer coefficients than laminar 
flow, any method of stabilizing the boundary layer is 
worth investigation. The statement by Lees! that 
ihere was a possibility of completely stabilizing the 
boundary layer was accordingly welcomed by aero- 
thermodynamicists. The numerical investigation of 
that possibility using the realistic Prandtl Number of 
0.75 and the Sutherland viscosity-temperature law is 
therefore the main purpose of the present paper. The 
Crocco method?* is employed to obtain the distribu- 
tions of flow properties across the laminar layer, and the 
Lees-Lin theory‘ is used in the stability analysis. The 
inviscid, as well as the viscous, solutions for complete 
stabilization of the laminar layer are obtained, includ- 
ing results for Prandtl Number unity and viscosity- 
proportional-to-temperature law. V-2 flight-test skin- 
temperature data are used to compare theory with ex- 
periment. Only two-dimensional flow on a flat plate 
with zero pressure gradient is considered. Prandtl 
Number and specific heat are assumed constant across 
the An interesting result is that, 


boundary layer. 
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Fic. 1. 
fluid. 


with zero pressure gradient, it is impossible to stabilize 


completely the boundary layer at speeds above approxi- 
mately Mach Number 9 with any amount of cooling. 


LAMINAR BOUNDARY-LAYER ANALYSIS 


Crocco’s Method for Flat Plates 

In order to carry out a stability analysis of a laminar 
boundary layer (Fig. 1), the flow-property distributions 
For a flat 
plate with zero pressure gradient and steady motion, 


within the layer must first be calculated. 


these distributions are obtained upon solution of the 
following respective equations of continuity, momen- 
tum, and energy for thin boundary layers: 


(0/Ox) (pu) + (0/dy) (pv) = 0 


Ou ‘ ou %) ( 2) (2 
= v>— = <) 
™ Ox ‘ Ov Ov 4 Oy 
rey) 4 rey) (=) 4 o (° *) (2) 
Uu v - = ° ”) 
: Ox F Ov “ Oy Oy \cp Oy 


along with the equation of state and a viscosity-tem 
perature relation. Owing to thinness, the pressure is 
constant across the layer, so that for a perfect gas the 
density varies inversely with temperature. Viscosity 
is assumed to be a function only of temperature. 

In these equations, u and v are the x and y components 
of the velocity at any point in the flow, the x-axis being 
taken along the plate in the direction of the free stream 
and the y-axis being taken perpendicular to the plate. 
The symbols p, u, k, p, ¢,, and 7 indicate the density, 
viscosity, thermal conductivity, pressure, specific heat 
at a constant pressure, and enthalpy per unit mass, 
respectively. 

Assuming that the Prandtl Number (Pr = c,u/k) was 
constant and that the enthalpy 7 was a function of ve- 
locity u alone, and upon elimination of v through the 
continuity equation, Crocco transformed Eqs. (1)—(3) to 
a pair of simultaneous ordinary differential equations 


in gx (us) and 7%(u%)—namely, 


£,8, + 2u,p,u, = 0 


(6,” + Pr(u.?/t..)lg, + (1 — Pr) 4,’g,’ = 0 


Schematic sketch of boundary layer in compressible 


and the primes refer to differentiation with respect to 
u,. Subscript o refers to the free-stream conditions. 
u(ou/Oy). 
Following Crocco, Eq. (5) is first integrated to yield, 
in terms of 7,’(0) and g,(0), 
' ) x 


7 2S aie u 
1,’(u,) = io)| & | ~ Pr ( - 
g,(u), 7 
g,(0) ]' " “mT go (0) VP"~! 
du, 6 
g,(u,) 0 g,\u), 


sais si ul 
1,(0) + 2,'(0)-I(u,) — Pr- 


and shear 7 = 


1,(u,) = -S(Uu, 7 


1 


” tx a (u) Pi l 
| | f ’ | du, 
0 g, (0) 
” x g (u) Pr l bt o (u.) 1-—P 
Ft.) = I k “| dus, f * 2 
0 g.(0) 0 g,(O) 


Since 7,(1) = 1, it follows from Eq. (7) that 


l : > 
1 — 2,(0) + Pr-J(1) -— (8 
I(1) 1 


substitution of which in Eq. (7) gives 


where 


I(u,) = 


ds 


1/(0) = 


; : I (u 
i.(O) — [7,(0) — 1] 
I(1) 


Uu [A say 7 Jia) 9 
1 I(1) 


When the specific heat is constant, then 


i(u,) = 


Pr ‘ 


u.7/4. = (4 1)M 


for1,, = c,7°,, so that Eq. (9 


and7, = 7, = T/T 


becomes 


I(u,) 
ra = 7, (0) — [7,(0) — 1] : 
(4) ¢ ) | T(1) 
Pr-(y — 1)-M | ee 7A) — Ju, | (10 
I(1) 


in which 7 is absolute temperature, y is the ratio of 
the specific heats at constant pressure and constant 
volume, and \/, is the free-stream Mach Number 
i.e., the ratio of the velocity of the free stream to the 
velocity of sound in the free stream. In the special 
case where Pr = 1, Eq. (10) reduces to the more widely 
known Crocco temperature-velocity relation 


T.(0) — {[7,(0) — 1]u, + 


? xa Yartu(t —u,) (ll 
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STABILITY OF THE 
The symbol 7,(0) will be frequently referred to as 
Tw/T., in the following discussion. 

In principle, Eqs. (4) and (9) can be solved for 7, and 
g, by successive approximation for arbitrary Prandtl 
Number once a viscosity-temperature relation is as- 
sumed. However, Crocco has shown, through ex- 
tensive numerical calculations, that for boundary- 
laver flow the enthalpy-velocity distribution is, for all 
practical purposes, independent of the viscosity-tem- 
perature law when the Prandtl Number is not too far 
fom unity. Hence, Eq. (9) represents immediately, to 
a good approximation, the enthalpy distribution upon 
substitution of a known shear distribution correspond- 
ing to any viscosity law. Accordingly, in order to sim- 
plify calculations, the viscosity-temperature law may be 
chosen such that Eq. (8) reduces to 


She” + 2u, = 0 (12) 


which is the Blasius equation for incompressible flow the 
solution of which is already tabulated.” * Such a 
simplificatiop requires that p,u, = 1. The /’s and J’s 
corresponding to the Blasius solution for g, are given 
in reference 2 and reproduced in reference 3. The 
enthalpy distribution can now be readily computed for 
arbitrary Prandtl Number, Mach Number, and wall-to- 
free-stream temperature ratio. Once the enthalpy 
distribution is calculated, an accurate solution for g, 
isnext obtained upon integration of Eq. (4) for the par- 
ticular viscosity-temperature law desired. 

Crocco, furthermore, developed an accurate iterative 
procedure for solution of the momentum equation 
[Eq. (4)]. Numerical analysis for g,(u,) using this 
procedure has been carried out in reference 3 for a wide 
range of Mach Numbers and wall-to-free-stream tem- 
perature ratios as a background for the present sta- 
bility analysis. A Prandtl Number of 0.75 was assumed, 
and the Sutherland viscosity-temperature law was 
used. From the data of Tribus and Boelter,® the Suth- 
erland S constant was found to be 198°R. for air in 
order to fit the viscosity law to the data starting with 
an initial ambient temperature 7, of 392°R., whence 
2/2... = 0.505. 

In recapitulation, the Crocco method provides a 
rapid means of calculating the enthalpy (temperature, 
density) distribution and a convenient numerical in- 
tegration process for solving the momentum equation 
to determine the shear distribution. With these two 
distributions, all other properties of the laminar bound- 
ary layer can be computed. 

In the calculation of the geometric distributions of 
the flow properties, the normal distance parameter used 


in this report is 


vy lo ux La 
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since r = y(du/dy). 


LAMINAR 


C 
+ 
Os 


BOUNDARY LAYER 


Cone Solution 


In the calculation of the stability of the boundary 
layer on cones, it is necessary to take into account the 
geometric effect upon the thickness of the layer. 
Hantzsche and Wendt® have shown that the equations 
for a thin laminar boundary layer on a circular cone 
at zero angle of attack in a supersonic stream with 
attached shock wave can be reduced to equations hav- 
ing the same form as those for a flat plate when the 
normal distance parameter (y/x) V Re, for the cone 
is taken as 1/+/3 times that for the flat plate. It fol- 
lows that, for the same local flow conditions just outside 
the boundary layer, the cone Reynolds Number cor- 
responds to one-third of its value on a flat plate. This 
result is valid for stability, as well as for coefficients of 
heat transfer and skin friction. 


STABILITY ANALYSIS 


The boundary-layer stability analysis is an eigen- 
value problem in that the relation between certain 
parameters is investigated which allows the boundary 
conditions to be satisfied by the solution of the under- 
lying differential equations. The parameters involved 
in the flat-plate laminar layer are the wave length \ of 
an infinitesimal disturbance, the wave velocity c, of the 
disturbance, the free-stream Reynolds Number Re,, 
the free-stream Mach Number J/,, the wall-to-free- 


stream temperature ratio 7/7, and the Prandtl 
Number Pr; thus, 
G = hh; Re, MU. Taft st) (13) 


The boundary conditions are on the velocity disturb- 
ances at the wall and free stream. The underlying 
differential equations are five linearized disturbance 
equations obtained from the five basic instantaneous 
equations for the two-dimensional motion of a perfect 
gas (two equations of motion, one equation of conti- 
nuity, one equation of energy, one equation of state) as 
follows: First, any quantity Q subject to fluctuation 
is separated into a steady-state part Q and a small dis- 
turbance Q’; thus, Q = 0+ Q’. After substitution of 
these sums into the basic equations, the mean tem- 
poral equations are subtracted out, quadratic terms in 
the small disturbances are then neglected, and the thin 
boundary layer assumptions (viz.,v <<< u, 00/ax«x « 
0Q/dy) are applied. The resulting equations are too 
lengthy to be repeated here, see reference 4. 

The disturbances will, in general, have a wavelike 
form and therefore can be separated by Fourier analy- 
sis into sine waves. Since the disturbance equations 
are linear in Q’, each sine wave can be treated sepa- 
rately. Thus a solution of the following type is investi- 
gated: 


ia(x—cl) 


(14) 


Q'(x, y, t) = glyje 


where q(y) is the amplitude of the disturbance and will, 
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incompressible flow (Blasius solution). (c) Neutral stability 


curve for incompressible flow. 


in general, be complex; a is the wave number defined 
by 2 w/\ and taken as real and positive; and c may be 
complex. The real part c, of c is the wave velocity, 
whereas the imaginary part c; indicates whether the 
oscillation is self-excited, neutral, or damped, accord- 
ing to whether the imaginary part is positive, zero, or 
negative. The problem of stability is therefore the 
determination of the sign of c; such that the disturbance 
equations and the boundary conditions are satisfied. 
In particular, in the present paper the minimum 
Re, for c, = 0 and arbitrary J, and 7yw/T,, is de- 
sired. For example, Fig. 2c represents the neutral 
stability curve for MW, = 0 and 7w/T, = 1 
Cooling will shift it to the right; 
Increase in Mach Number when 


(incom- 
pressible flow). 
heating, to the left. 
the plate is insulated will also shift it to the left. 


Inviscid Solution 


Valuable information on boundary-layer stability 
can be obtained by first completely neglecting the 
effect of viscosity in the linearized disturbance equa- 
tions. This gives the inviscid solution. Thus, in this 
type of solution only pressure and inertial forces are 
assumed acting on the disturbances in the established 
(by viscosity) shear flow. As shown by Lees and Lin, 
the inviscid case is actually the limiting case of infinite 
Reynolds Numbers. 

For boundary-layer profiles in incompressible fluids, 
Rayleigh and Tollmien concluded, for the inviscid case, 
that the necessary and sufficient conditions for the 
existence of neutral and self-excited disturbances are 
that the velocity profile has an inflection point between 
the plate and the free stream. They also showed that, 
in the case of a neutral disturbance, c = u at the in- 
flection point. Since the Blasius profile has no in- 
flection point (see Figs. 2a and 2b) although d*u/dy? = 
0 at y = O, such profiles are always neutrally stable 
regardless of Reynolds Number. The neutral curve of 
Fig. 2c then moves to the right to infinity or, rather, 
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does not exist for the inviscid case with zero pressure 
gradient. An adverse pressure gradient would cause an 
inflection point to appear in the profile, and therefore 
the boundary layer would be immediately unstable for 
any Reynolds Number so that the neutral curve of Fig, 
2c would then be the vertical ordinate through Re = 9, 
A favorable pressure gradient is, on the other hand, 
stabilizing. Although these results are valid only for 
the limiting case Re > ~, they do offer some informa- 
tion in that they indicate the destabilizing effect of a 
point of inflection in the velocity profile of a real in- 
compressible fluid. 

Lees and Lin have extended the Rayleigh-Tollmien 
criteria for the existence of neutral and self-excited dis- 
turbances to compressible fluids. They have shown 
for boundary-layer profiles that the necessary and 
sufficient conditions for the existence of neutral subsonic 
disturbances is that the quantity (d/dy)[p(du/dy)| 
vanishes for some value 


“u.>1— (1/M,) 


Ue = = C 


They have also shown that this is a sufficient condition 
for the existence of self-excited subsonic disturbances, 
Disturbances that are subsonic relative to the free 


stream—i.e., 


ce ] 
t ~—-¢<@. 2 >1- 
u. M 


x 


where a is sound velocity—are called subsonic dis- 


turbances and are shown to be the only kind having sig- 
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nificance in boundary-layer stability because they exist 
only for discrete eigenvalues of c, a, Re, M.., Tw/T.., 
and Pr. Sincec, = u, at the point where 


(A) 


it follows that in subsonic flow a neutral or self-excited 
subsonic disturbance exists whenever Eq. (A) holds 
at any point in the layer. It is to be noted that, in 
general, the distribution of p,(du,/dy,) is decisive in 
the determination of the stability of the laminar layer. 
The quantity d/dys[p,(du,/dy),] reduces to d*u,/- 
dy,? when the fluid is incompressible such as in the case 
of liquids or with low-speed gas flow without heat 
transfer. However, with heat transfer in gases, the 
quantity (d/dy,) [p,(du,/dy,)] must be considered even 
though 7, — 0 (see Fig. 3). 


(d/dy,) [p,(du,/dy,)] = 0 


According to Lees and Lin, the point at which Eq. 
(A) holds for u, = 1 — (1/M,) divides stability from 
instability as controlled by Mach Number, wall-to-free- 
stream temperature ratio, or Prandtl Number. This 
point is important because of its simplicity of calcula- 
tion and its approximation to the more exact conditions 
when viscosity is included. 


The process by which cooling stabilizes the boundary 
layer is shown in Figs. 3 through 6. For M, = 0, it 
is seen that the layer is unstable [bulge in p,(du,/dy,) 
away from the wall] according to inviscid theory when- 
ever the layer is heated; vice versa, the layer is stable 
(no bulge) whenever the layer is cooled. At M, = 4 
and 8, the temperature ratio required to stabilize the 
layer completely is found at about 2.4 and 1, respec- 
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tively, as indicated by the crossing of the hump in 
p,(du,/dy,) at the point u, = 1 — (1/M,). When 
M., = 12, it is seen that the boundary layer can never 
be completely stabilized by cooling regardless of Reyn- 
olds Number because the point at which Eq. (A) holds 
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Fic. 5. Distribution of product of density and vorticity across 
laminar boundary layer for ., = 8 and various wall-to-free- 
stream temperature ratios. 
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Furthérmore, dux/dyx = dT %/dyx since we = Tx. 
Therefore, the calculations are considerably simplified 
when peux = 1. 
(3) Assumptions: 
proportional to temperature; specific heat constant. 
For Pr | and cy = constant, the simple relation 


between temperature and velocity given by Eq. (11) is 


Prandtl Number unity; viscosity 


used, and the brackets involving gs yield unity in 
Eq. (17). Furthermore, since pxux 1, the shear dis 
tribution is the Blasius distribution. Therefore, the 
calculations for this case reduce to a simple set that 
can be carried out with as much accuracy as desired. 

In Table 1 are listed the results of the three sets of 


calculations for the inviscid solution. 


TABLE | 
Temperature Ratios Required for Complete Stabilization of the 
Laminary Boundary Layer Based on the Inviscid Solution 


§ = Py Pr = Pr 
0.505 0.75 Pate = 1 0.75 Pe bn l 1.0 
Tw/1 VW Tw/T VW Tw/1 VW 
l 0) l 0 l 0 
1.17 l 1.18 l 1.17 l 
1.38 1.5 1.35 1.5 1.25 1.25 
1.65 2 1.565 2 1.34 1.5 
1.91 2.5 1.74 2.5 1.438 1.75 
2.14 5 1.86 3 1.50 9 
2.41 ! 1.93 } 1.55 2.5 
40 5 1.74 5 1.55 3 
» 2) 6 lL. 2 6 1.31 } 
1.75 7 0.45 7 0.72 5 
1.08 8 0) 7.48 0 5.8 
() 8.97 


Viscous Solution 

The viscosity terms are now considered to a first ap- 
proximation—1.e., Reynolds Numbers are large but 
finite. Just as Lin’ had modified Heisenberg’s criterion 
for viscous incompressible fluids to read, “every sym- 
metrical or boundary-layer profile is unstable for suffi- 
cently large values of the Reynolds Number,” so 
has Lees' extended Heisenberg’s criterion to compres- 
sible fluids for subsonic disturbances. Lees concluded 
that when 7. < 1 any laminar boundary-layer flow in 
a viscous conductive gas is unstable at sufficiently high 
(but finite) Reynolds Number. For M, > 1, Lees 
further showed that any laminar layer for which the 
quantity (d/dy) [p(du/dy)| vanishes for some value of 
ux > 1 — (1/M_.) is also unstable at sufficiently high 
but finite) Reynolds Number. When (d/dy) [p(du/- 
dy)] is negative for ux > 1 — (1/M.,), the laminar layer 
will still be unstable at sufficiently high Reynolds Num- 
ber until (d/dy) [{p(du/dy)] becomes large enough neg- 
atively as controlled by the free-stream Mach Number 
and the rate of heat transfer from the boundary layer 
to the wall. Thus, the effect of viscosity is destabilizing 
at high Reynolds Numbers. On the other hand, since 
the flows will become stable for sufficiently low Reynolds 
Numbers, the neutral stability curves for incompres- 
sible, as well as compressible, flow are expected to be of 


the type indicated in Fig. 2c. Hence, a whole range of 
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minimum critical Reynolds Numbers will exist on either 
side of the inviscid solution curves of Fig. 7. When 
M.. < 1, the boundary layer will be unstable for suffi- 
ciently high Reynolds Numbers regardless of cooling 
rate. However, when 1/. > 1, the rate of cooling can 
be great enough to stabilize the boundary layer com- 
pletely regardless of Reynolds Number, at least over a 
range of Mach Number, in which case the curve of Fig. 
2c is shifted again to the right to infinity. 

The procedure used in this report for calculating the 


neutral stability curve a = a(Re, M., Tw/T., Pr) is 
that given by Lees and Lin in reference 4. It was 


shown in reference 4 that for the viscous case, the rela 
tion between the parameters takes the form 
(1 + A) (U + 1V) 
?(Z) t : : - (23) 
1+ A(U 4+ 1V) 
where ® (Z) = [1 — F(Z)]~!, in which F(Z) is the Tiet 
jen’s function defined by 


F(Z) 1 + = 3 (24) 
Zz / ¢ H, Kop Ja 
Z = | =e (du el) a (95) 
Vx r 
Furthermore, A is defined by 
(dix dy JwV* cy (26) 


A(cs) = 

Ox 
and U and V are real functions of ax, cx, 7.., Tw/T., 
and Pr. The quantity /7,;"" is the Hankel function of 
the first kind of order '/;. In Eq. (25), the quantities 
dux/dyx, vx, and yx are calculated at ux Cx. Also, 
2rL/X and Rex. pt. L/a., 
where the reference length 


by definition, ax = 


L = x/Vp.u.x/p 


This reference length L is used in a and Re because it 


is used to normalize y in 


2.6 0.8 ,— . . . ay r 
| 
2.2 os——. ‘ { } | 
| ©. (2) 
1.8 0.4 }— 4 4 4 | 
®, (2) 
1.4 0.2 { { } } 


@(2) (2) 
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Fic. 8. Real and imaginary parts of the function #(Z) = 1/ 
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yx = (y/x) V puX/M, 
It then follows that 


Re,., — PU aX He - (Rex ..)” 


in which x is the distance from the leading edge of the 
plate. 

The real and imaginary parts of $(Z) have been com- 
puted in reference 4 and are reproduced in Fig. 8. 

Lees and Lin give as an approximation for V the 
following formula, which has been justified numeri- 
cally”! for incompressible flow and low supersonic 


speeds: 
(**) 
- 
; dysiw [| Te? a ( og 
as = : ao px 
7 *W a)" dyx dy 
dy x lg = Ce 
(27) 


in which the subscript W refers to the wall. 
Separation of the real and imaginary parts of Eq. (23) 
yields 
(Z) = (1+ A) [U(1 + AU) + AV?] X 
((1 + AU)? + A?V?] 


®(Z) = (1 + A)V[(1 + AU)? + A?V?]- 


(28) 
(29) 


Since A and V are given by Eqs. (26) and (27), U may 
then be determined from Eqs. (28) and (29) such that 
®,(Z) and ®,(Z) are satisfied at the same value of Z in 
Fig. 8. The procedure is rapid. <A and V are first 
plotted from Eqs. (26) and (27) as curves versus, say, 
Tw/T.. Sliding along these A and V curves, U can 
be successively computed from Eq. (29), for a chosen 
@,(Z), and substituted into Eq. (28) until the corre- 
sponding value of ®,(Z) is found. 

In order to obtain a as a function of Rex ., the pro- 
cedure is first to calculate the product axRex.,.. For a 
given set of 1/,, 7w/T., and Pr, the quantity ax Rex , 
can be determined from Eq. (25) as a function of Z. 
Values of ux = cx [greater than 1 — (1//,)] at which 
Eq. (25) is evaluated are, however, first determined 
such that the corresponding values of A and V from Eqs. 
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(26) and (27) satisfy Fig. 8 at each selected Z. The 
calculation of ax at the corresponding value of Z would 


then give the eigenvalue for Rex. from ax Rey... 


The determination of the minimum critical (neutra] 
Reynolds Number would require a careful calculation 
of @ as a function of Z. However, the calculation of 
a is not easy and was not undertaken in this investiga. 
Rather, the criterion of Lin, 
that the minimum critical Reynolds Number corre. 


tion for lack of time. 


sponds to the peak of the ®,(Z) curve, and an empirical 


formula of Lees [Eq. (30) ] were used. 





That the minimum critical Reynolds Number cor. 
responds well with the peak of the ®,(Z) curve for 
incompressible flow is seen from Fig. 9 in which are 
plotted detailed calculations carried out by Lin.’ In 
this figure the subscript 1 indicates that the reference 
length used was the displacement thickness 6*; also, 
all upper branches correspond. The upper branch of 
the ®;(Z) curve in Fig. 9 corresponds to the right-hand 
branch of @,(Z) in Fig. 8. The peak of the $,(Z) curve 
is 0.58, and the corresponding ®,(Z) is 1.48 at Z = 3.22. 
It is assumed in the present paper that this criterion 
also holds for compressible flow. 

Since it is shown for supersonic flow! that ax — 0 as 
Cx > 1 — (1/M,), it follows that Rex . 
Z or, rather, finite ax Rex .. in Eq. (25). 
the possibility of finding a wall-to-free-stream temper- 


— o for finite 
Hence, there is 


ature ratio for which the boundary layer is completely 
stabilized regardless of Reynolds Number in supersonic 
flow. Thus, for each Mach Number /_, the tem- 
perature ratio 7j/7, was desired at which ®,(Z) = 
0.58 and ®,(Z) = 1.48 at ue = 1 — (1/M,) using 
Eqs. (26) to (29). The results of such calculations are 
shown in Fig. 7 for the three sets of conditions used in 
the inviscid solution. The numerical results for the 
viscous case are listed in Table 2. 

It will be noted that the curve for the case Pr = | 
and pxusx = | is different from that previously given 
by Lees! who assumed A in Eqs. (28) and (29) to be 
small; however, this is not the case for the higher Mach 
Numbers and lower temperature ratios (see Fig. 15). 


TABLE 2 
Temperature Ratios Required for Complete Stabilization of the 
Laminar Boundary Layer Based on the Viscous Solution 


= Pr = Pr = Pr = 
0.505 0.75 Pele = 1 0.75 Pete = 1 1.0 
Tw/T« M. Tw/T. WU Tw/T M 
0 1 0 1 Q ] 
0.80 1.25 0.80 1.26 0.84 1.25 
0.99 1.5 0.99 1.5 1.01 1.5 
1.27 2 1.24 2 1.20 2 
1.51 2.5 1.43 2.5 1.27 2.5 
1.71 3 1.55 3 1.29 3 
1.95 4 1.65 4 L.a3 4 
1.99 5 1.53 5 0.64 5 
1.88 6 1.10 6 0 5.8 
1.54  § 0.41 ‘4 

0.93 8 0 7.48 

0 8.97 
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Fic. 10. Effect of Prandtl Number and viscosity-temperature 
law on the cooling required for complete stabilization of the lami- 
nar boundary layer. 


Fig. 7 indicates that, for the realistic flight conditions 
of Pr = 0.75 and Sutherland viscosity law (6 = 0.505), 
it is impossible to stabilize the laminar layer, com- 
pletely, above approximately M., = 9 with any amount 
of cooling when the pressure gradient is zero. That 
this is expected at high Mach Numbers can be observed 
in Fig. 6, where the hump in the ps(du«/dys) curve re- 
mains outside of the point uw. = 1 — (1/M,). [Since 
Vis positive from Eq. (29), it follows from Eq. (27) that 
(d/dy) [p(du/dy)| is always negative—i.e., corre- 
sponds to a point on the right side of the hump in 





p(du/dy).] 

The relative effects of the viscosity-temperature law 
and Prandtl Number on the cooling required for com- 
plete stabilization of the boundary layer are shown in 
Fig. lO for Tw/T,, = 0. It is seen in the figure that at 
Pr = | the effect of variation in the viscosity law is not 
(Additional calculations for the entire viscous 
1 and Sutherland law showed a rise to 
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about 7/7, = 1.5 compared with about 1.3 for the 
case of Pr = 1 and peux = 1.) For Pr = 0.75, the 
effect of viscosity law is considerably greater than for 
Ir = 1. 

A formula given by Lees! was used to calculate the 


minimum critical Reynolds Numbers. The formula is 


25 (Tux), dux 
(30) 
cet» V1 — M2 (1 — ce)? \d¥e/ 


Re... = 
which fits results carried out accurately by Lees for 
M.. = 0.7 and variation of wall-to-free-stream tem- 
perature ratio from 0.70 to 1.25 using the power vis- 
cosity law with exponent equal to 0.76. Although Eq. 
(30) can hardly be expected to hold for large variations 
in Mach Number and temperature ratio, yet it indi- 
cates the trend of the minimum critical Reynolds Num- 
ber and therefore was used to obtain Fig. 11. How the 
minimum critical Reynolds Number varies with tem- 
perature ratio is shown in Fig. 12 for M, = 4. It is 
immediately apparent from Figs. 11 and 12 that the 
cooling rate itself is critical near the infinite Reynolds 
Number curve, because a small change in temperature 
ratio produces a large change in the neutral-stability 
Reynolds Number. This result is important in a prac- 
tical sense, since it indicates a radical change in the 
state of the boundary layer as the curve of complete 
stability is approached. Fig. 11 also shows that, above 
about 7, = 9, the laminar boundary layer with a zero 
pressure gradient is in a state of impending transition 
to turbulence. 

In order to give some idea of the effect on the sta- 
bility of the laminar boundary layer of cooling due to 
heat radiation from the surface, the temperature of a 
plate with emissivity unity was calculated and the 
wall-to-free-stream temperature ratio was plotted in 
Fig. 11. It is evident that even full radiation cannot 
completely stabilize the boundary layer, contrary to 
what had been predicted earlier.' 
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Minimum critical Reynolds Number as a function of free-stream Mach Number and wall-to-free-stream temperature 
i Prandtl Number 0.75 and Sutherland viscosity law. 
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STABILITY OF TAs 


Fig. 13 gives the dimensionless wave velocity c¥ as a 
‘ynction of Mach wall-to 
jree-stream temperature ratio, all for Pr 0.75 and 
0.505), such that &,(Z) 


free-stream Number and 
Sutherland viscosity law (@ 
(58. 

For future work, 1n case ay is calculated, the product 
vx Rex 
ind was plotted in Fig. 14. 


299 


was computed from Eq. (25) for Z 3. 


Fig. 15 shows the magnitude of A of Eq. (26) for a 
that A cannot always be 


neglected in compressible flow, although it may be 


few cases. It is evident 


neglected as was justifiably assumed by Lin for incom 

pressible flow. 
EXPERIMENTAL VERIFICATION-USE OF STABILITY 

RESULTS IN AERODYNAMIC HEATING ANALYSIS 


It has been predicted that heating destabilizes and 
cooling stabilizes the laminar boundary layer (Fig. 12). 
\Ithough this conclusion was based on an analysis as 
suming infinitesimal disturbances, nevertheless it is ex 
pected to apply generally to the advancement or delay 
That heat 


ing advances and cooling retards transition has been 
13 


of transition of laminar to turbulent flow. 


verified experimentally by a number of investigators.® 
Fig. 16 shows some data from Higgins and Pappas!* ob 
tained by observation of surface-tube Mach Number 
readings as a function of Reynolds Number through the 
transition region. The data indicate a marked influ 
ence of heating on decreasing the Reynolds Number of 
transition. 

Flight-test data have also corroborated the effect of 
heat transfer on the state of flow of the boundary layer. 
Of particular interest are the data reported by Fischer 
and Norris'* on surface-temperature measurements ob 
tained on the smooth conical nose of a V-2 rocket dur 
ing ascent in New Mexico on October 9, 1947. Fig. 17 
shows the skin station G 
located 12 in. from the tip of the nose cone. 


temperature measured at 
There was 
no boundary-layer trip before this station. For the 
purpose of comparing the data with theory, theoretical 
temperature curves for laminar and turbulent layers 
ire also plotted in Fig. 17. It is evident from the 
figure that transient cooling of the boundary layer has 
stabilized the flow, as signified by the switch of the data 
from the turbulent-flow to the laminar-flow theory. 
The course of events is best seen when the data are 
plotted on a minimum critical Reynolds Number dia- 
gram as in Fig. 18. (Since Fig. 18 is a portion of Fig. 
\l for flat plates, one must bear in mind that the Reyn- 
olds Numbers indicated on the diagram are equivalent 
to three times their values when applied to cones, in 
accordance with the previous remarks about the rela- 
tion between laminar boundary layers on cones and 
flat plates.) Also shown in Fig. 18 is the temperature 
required to stabilize completely the boundary layer at 
the measuring station. This latter curve is obtained 
by plotting one-third of the actual local cone Reynolds 
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Number at the station in question. Where the two 
temperature curves (actual and required) intersect 
corresponds to the instant of complete stabilization up 
to that point of the missile. The instant of complete 
stabilization is seen to occur at a local cone Mach 
Number of 2.05, and the corresponding time is shown 
in Fig. 17 where it appears to indicate the beginning of 
the strong deviation toward the laminar-flow theory. 
Thus the fundamental theory of stability of small sub- 
sonic neutral disturbances in compressible flow boundary 
layers seems to have some substantiating evidence in 


its favor. 
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7 Mach Number depending upon Prandtl Number anq 
viscosity-temperature law. 
ACTUAL > (2) For Prandtl Number of 0.75 and Sutherland yijs. 
LOCAL Re nn cosity law, the laminar nano epa om with zero pres. 
STABILIZATION sure gradient is unstable for M, > 9 with any amount 
6 of cooling. 
(3) The effects of Prandtl Number and viscosity. 
8 temperature law on stability are appreciable. 
e (4) Wind-tunnel experiments have verified qualita. 
© tively the effect of heating and cooling on stability. 
S (5) V-2 flight-temperature data have indicated sta- 
. bilization due to transient cooling. 
MINIMUM REFERENCES 
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1 2 lation, AL-684, North American Aviation, Inc., Los Angeles, 
LOCAL MACH NUMBER,M. Calif., 1948. 
Fic. 19. Graph for use in aerodynamic heating analysis. 3 Van Driest, E. R., Investigation of Laminar Boundary pr 


A better way of observing the course of events as far 
as stability is concerned is to plot the minimum critical 
Reynolds Number corresponding to the actual skin 
temperature (in the present example, using the above- 
mentioned V-2 data, three times the values obtained 
from Fig. 18) and the actual Reynolds Number cor- 
responding to the station in question against local free- 
stream Mach Number, as shown in Fig. 19. From such 
a diagram the advancement or retardation of turbu- 
lence can supposedly be seen as the curves diverge or 
converge with increase in Mach Number. The inter- 
section of the curves again corresponds to the instant of 
complete stabilization. Comparing Fig. 19 with Fig. 
17, one would find that the region of convergence in 
Fig. 19 corresponds roughly with the interval of di- 
vergence of the data and the turbulence theory be- 
tween about 38 and 46 sec. In addition to its use in 
the interpretation of temperature data, such a plot as 
Fig. 18 has practical use in the calculation of skin tem- 
peratures of high-speed vehicles because it would indicate 
to the designer the change of state of the boundary layer 
as the vehicle moved along its trajectory. 


For comparison with data at station G, skin tempera- 
tures were also measured at station H on the opposite 
side of the nose cone during the above-mentioned V-2 
flight. However, a boundary-layer trip had been in- 
stalled upstream of that station. That the data fol- 
low the turbulence theory well when stabilization of 
the boundary layer by transient cooling is prevented 
by means of a trip is evident from Fig. 20. 


CONCLUSIONS 


(1) As predicted by Lees, the laminar boundary 
layer in supersonic flow with zero pressure gradient can 
be completely stabilized by cooling regardless of Reyn- 


olds Number, however, only over a limited range of 
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Theoretical Determination of the Minimum 
Drag of Airfoils at Supersonic Speeds 


ROBERT T. JONES* 


Ames Aeronautical Laboratory, N.A.C.A. 


ABSTRACT 


Consider a thin wing in frictionless flow and suppose the plan 
fam of the wing and also the total lift to be given. The drag of 
the wing will then depend on the way in which the lift is dis- 
tributed over its surface. Ina previous paper, it was shown that 
the minimum drag occurs when the superposition of the induced 
disturbance fields in forward and reversed motion results in a 
constant value of the induced downwash at all points of the wing 
surface. Similar problems involving the ideal distribution of 
thickness over the surface were found to lead to similar condi- 
tions governing the distribution of pressure in the superimposed 


or “combined”’ flow field. 

The present paper describes a method for determining mathe- 
matically the combined disturbance field, and in certain cases the 
minimum drag, of wings at supersonic speeds. ‘The simplest 
analytic example is provided by the wing of elliptic plan form, 
which achieves its minimum drag when the lift is distributed 
uniformly over the surface. With a symmetrical distribution of 
thickness, the requirement of minimum drag for a given total vol- 


ume is found to lead to profiles of constant curvature 


INTRODUCTION 


y THE THEORY OF WINGS at subsonic speeds, it is 
shown that the production of lift by a wing of finite 
span gives rise to a drag force that depends on the dis- 
tribution of lift over the span. This component of the 
drag, which arises in frictionless motion, may be re- 
lated to the energy required for the continual extension 
of the two-dimensional field of motion induced by the 
wake of trailing vortices. Alternatively, by examining 
conditions in the vicinity of the wing sections, the drag 
may be related to the downward inclination of the air 
stream induced at the position of the wing by the ac- 
tion of the trailing vortices. Following the latter con- 
cept, the drag arising from the lift at subsonic speeds 
has been termed the “induced drag.”’ It was shown by 
Munk' that this drag is a minimum for a given lift and 
a given span when the induced downwash is constant 
at all points of an equivalent lifting line, or vortex, hav- 
ing the same spanwise distribution of lift as the wing. 
It was further shown by Munk that the induced drag 
is actually independent of the chordwise distribution of 
lift. 

At supersonic speeds an additional component of 
drag arises because of the formation of waves by the 
airfoil, and in this case the drag depends on both the 
spanwise and chordwise distributions of lift or, in other 
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words, on the actual distribution of lifting pressure 
over the surface of the wing. In order to extend Munk’s 
problem to wings at supersonic speeds it was necessary 
therefore to consider not merely the span as given but 
the actual shape of the wing in plan view. The prob- 
lem could then be stated in the following form: For a 
given plan form s and a given total lift 1, what dis- 
tribution of the lift Z over the surface s results in the 
minimum drag? 

Furthermore, at supersonic speeds a certain drag 
arises from the thickness of the airfoil independently 
of the lift. The two components of drag may, how- 
ever, be considered separately and later added in any 
desired combination. To isolate the effect of lift, as dis- 
tinct from the effect of thickness, it is sufficient to re- 
place the wing by its mean surface, which is supposed 
to be warped or cambered in whatever way may be re- 
quired to cause the specified distribution of lift. On 
the other hand, the drag arising from the thickness may 
be determined by considering the thickness to be sym- 
metrically disposed above and below a flat mean sur- 
face having no lift. In this way additional problems 
involving the ideal distribution of thickness over the 
plan form become apparent. 

In reference 2 it was shown that all distributions of 
lift having the minimum drag for a given plan form and 
a given total lift are characterized by a single condi- 
tion. If we suppose the wing with its given distribution 
of lift to be held fixed in a stream of velocity V, then 
there will arise in the vicinity of the wing and its wake 
additional small disturbance velocities u, v, and w. 
Now let the direction of the stream be reversed, but sup- 
pose that the curvature and inclination of the surface 
is so modified as to maintain the original distribution of 
lift. A new field of disturbance velocities u, v, and w 
will appear. The wake of trailing vortices will have 
the same form as before, but the wake will now extend 
from the wing in the opposite direction. We may now 
superimpose the two fields of disturbance velocities 
and obtain by this means a ‘“‘combined disturbance 
field,’ associated with the given distribution of lift. 
For the drag to be a minimum, the downwash in the 
combined disturbance field must be constant at all 
points of the wing surface s. 

If the velocity of flight is subsonic, the superposition 
of the two fields can be shown to result in a two-dimen- 
sional field of motion identical in form to the velocity 
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Fic. 1. Lifting surfaces with superimposed disturbance fields. 


field of the vortex wake. For distributions of lift 
having the minimum drag, the downwash induced by 
the vortex wake in its own plane is a constant, and it 
will be evident that in this case the downwash is con- 
stant not only within the plan form of the wing but at 
all points of the vortex ribbon ahead of and behind the 
wing (Fig. 1A). 

At supersonic speeds the combination of the forward 
and reversed disturbance fields again produces an in- 
finite, parallel vortex ribbon, but the field is no longer 
two-dimensional in character and is bounded by two 
overlapping zones of influence or wave fronts, as illus- 
trated in Fig. 1B. 

By treating the problem of thickness in a similar 
manner, it was also shown in reference 2 that the mini- 
mum drag for a given frontal area of the wing occurs 
when the pressure in the combined flow field is constant 
at all points of the wing surface. Similarly, consider- 
ation of the minimum drag consistent with a given total 
volume led to the requirement of a constant streamwise 
gradient of the pressure in the combined flow. 

The present paper describes a method for determin- 
ing the combined disturbance fields associated with 
given distributions of lift or thickness. The basic idea 
of the method is to represent the elementary solutions 
of the flow equation, such as the solutions for the source 
and for the horseshoe vortex, by contour integrals, 
following forms introduced by Whittaker*® and Berg- 
man‘ rather than the usual forms. The distribution 
of sources or dipoles over the wing surface is then rep- 
resented by a triple integral, in which the surface in- 
tegral, after a change in the order of integration, repre- 
resents a distribution of two-dimensional disturbances 
over the surface. The three-dimensional flow is thus 
obtained finally by the superposition of elementary 
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two-dimensional flows. As will be shown, this method 
enables the calculation of three-dimensional wing flows 
that satisfy the conditions for minimum drag and pro. 
vides, as examples, formulas for the minimum drag 
of wings of elliptic plan form at supersonic speed. 


PRELIMINARY CONSIDERATIONS 


As is well known in the thin-airfoil theory, the lift 
distribution over a thin cambered wing or lifting sur- 
face appears as the resultant of two equal, opposite 
pressures over the upper and lower surfaces. With 
the pressure disturbance given by 


Ap = —puV (1) 


where u is the longitudinal perturbation velocity, we 
obtain, for the local lift, 


l(x,y) = 2pVu(x,y,z); 2—~> +0 (2 


The velocity u is discontinuous across the lifting sur- 
face and is to be evaluated on the upper side (z > +0). 
In the lifting case, the downwash velocity w(x,y) is 
continuous over the whole plane of the wing. The 
drag is given by 


P w(x,y) 
p= f f l(x,y) Vv dx dy (3) 


In certain cases, the lift density / and the stream in- 
clination w/V may approach infinite values around 
the edges of the surface. The integral (3) must then be 
evaluated by a suitable limiting process, as described 
in reference 5. 

It has been shown by von Karman® and Hayes’ 
that the drag of a given distribution of lift is unchanged 
by a reversal of the direction of motion. Hence, the 
drag of a specified distribution of lift may be calculated 
from the corresponding distribution of downwash in 
either direction of motion or from the combined down- 


wash @, as indicated by the following formulas: 


I a l _ 
Dp = 7 ij lw dx dy = 7 ¥i lw dx dy = 
— 
: lf lw dx dy (A) 
2V Js . 


For the minimum drag the combined downwash @ will 
be constant over s, and we have 


D = (1/2)L(@/V) (5) 


where L is the total lift. 

Turning to the case of a prescribed distribution of 
thickness with no lift, it is noted that the velocity u and 
the pressure as given by Eq. (1) are continuous across 
the upper and lower sides of the mean plane, but the 
velocity w has a discontinuity related to the equal and 
opposite slopes of the upper and lower wing surfaces. 
If ¢(x,y) is the prescribed thickness of the wing at the 
point (x,y), we have 








wh 
the 
pot 
hov 


dra 
nes 
fort 


or 


ove 


air 
fyi1 


Th 
fiel 
bin 
Ma 
M 


sult 


wh 


a, | 


of t 


(Se 


\ 


she 


lethod 
flows 
d pro- 
| drag 


he lift 
g sur- 
posite 

With 


(1) 


Y, we 


r sur- 
+(), 
VY) is 

The 


(3) 


n in- 
ound 
-n be 
ribed 


ryes’ 
nged 
, the 
ated 
h in 
ywn- 


1 of 
and 
ross 
the 
and 
ces. 
the 











MINIMUM DRAG OF 


w(x,y) = —(1/2)Vt'(x,y); 2—~> +0 (6) 


where ¢’ denotes dt/dx. It can be readily verified that 
the velocity @, in the combined flow field, vanishes at all 
points of the wing plan form. A value of # remains, 


however, and determines the drag through the relation 


vec 
D=- *. / fri dx dy (7) 


As shown in reference 2 the requirement of minimum 
drag with various specifications on the maximum thick- 
ness or the volume of the wing leads to conditions of the 


form 


= constant (8) 


>, 


or 
Ou/Ox = constant (9) 


over all or part of the wing plan form. 


GENERAL FORM OF THE SOLUTION 


The field of disturbance velocities surrounding the 
airfoil will be characterized by a velocity potential satis- 
fying the well-known differential equation 


¢:: = 0 (10) 


(iM? — 1) Grr — Py ~ 


The same differential equation holds for the disturbance 
field in either direction of motion, as well as for the com- 
bination of the two fields. Since the fields for different 
Mach Numbers differ only by an affine transformation, 
it will be convenient to perform the calculations for 
M = +/2. 

As may be shown by direct differentiation, the re- 
sulting equation possesses the primary solutions*® 


¢ = F(ax — By — 72) (11) 


where F is an arbitrary, differentiable function and 
a, 8, and y are parameters determined so that 
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e&—-BP-—-y7y=0 (12) 


Through Eq. (12), a, 8, and y may be made to de- 


pend on a single complex parameter A. Writing 


\ = e” and setting 

a=] | 

B cos 6 = (1/2) [(1/A) + A] (13) 
vy = sin 60 = (1/2) [(1/A) ~ y)) 


Eq. (12) is satisfied for values of \ extending over the 


entire complex plane. For real values of @ (|A| = 1), 
Eq. (11) becomes 
yg = F(x — ycos@ — zg sin 6) (14) 


and the solution is seen to represent a plane wave of 
arbitrary form F. The wave front lies at an angle of 
45° to the x axis (MZ = +/2) but is inclined at an angle 
6 in the y,z plane. On the other hand, for large values 


of \ we have 
g = F[-—(A/2)(y — 22)] (15) 


and the solution here represents a cylindrical or two- 
dimensional flow with its axis parallel to x. This two- 


dimensional field is evidently a solution of 
“he te = 0 (16) 


with ¢,, separately equal to zero. 

A general solution of Eq. (10) may be constructed by 
superimposing a number of solutions of the form (11) 
for various values of the parameter \. It is clear that 
the form of the function F need not be the same for all 
values of \, so that F may depend on the two variables 
ax — By — yzandX. Thus we obtain 


g = f F(ax — By — yz,A) dr (17) 


where C is some contour in the \ plane. Eq. (17) is 
closely analogous to Whittaker’s solution of the La- 
places equation*® and belongs to the more general class 


oi integral operators studied by Bergman.‘ 


EXPRESSIONS FOR ELEMENTARY DISTURBANCE FIELDS 


As an example of Eq. (17), we may construct the well-known solution for the supersonic point source by means 


of the superposition of plane waves. 


This solution can be represented in terms of real values of 6 as follows: 


0, unless x? > y? + 2? 


forx > 0 


1 2 dé — 18 
D. = ‘ - “ “ 
7 ir? Jo x — ycosé — zsin#é si | 


(See references 9 and 10.) Here, R = Vix? 


R = 0 represents the Mach cone, which extends both ahead of and behind the point source. 


y? — 2° and is assumed to have a positive real part. 





—1 
forx < 0 
2rR 

The equation 
The integral (18) 


shows a ‘‘zone of silence”’ in the space between the fore cone and the rear cone. 



































: z with the real part of R positive. The integral (19) may 
now be written 
a P 4 f ar Fe ad 
2 a Cc B if A adil 271? (y — 12) (A — e1) (A — 2) \<9) i. 
™ a To evaluate the integral (23) by the method of regj. gap | 
mh, Ps dues, it is necessary to investigate the positions of the rearW 
4 x poles «, and ein the plane as functions of the coordi- the o1 
MACH ie My nates x,y,z. Fig. 2 illustrates this correspondence, secon 
CONE \Y/ a The three-dimensional x,y,z manifold is represented on at in 
a \ the two-dimensional complex plane by the identifica the c’ 
o “ tion of points with rays. This representation was used Each 
by Busemann in his conical-flow theory.'! Each ray then 
XYZ SPACE drawn from the origin in (x,y,z) space appears as two ray 
points ¢, and ¢.in the \ plane. For rays drawn toward The 1 
A the positive direction of x and lying inside the Mach then 
cone (y? + 2? < x’), the point e«, will lie inside the unit 1/\A 
circle, |A! < 1, while the point e will lie outside this the c 
6 8 circle at the reciprocal radius—that is, and I 
a " Ie, (24 of th 
€2. 
z Investigation of Eq. (22) shows that for each point Ne 
o (x,y,z) in the space outside the Mach cone (y* + 2? > conte 
Cc 6 x”), both e, and €: lie exactly on the unit circle. Thus, the s 
Cc €, and ¢€ project the space (x,y,z) inside the Mach live | 
cone on a surface, while the whole portion of the space enclo 
D (x,y,z) in the ‘‘zone of silence’ outside the Mach cone 
is represented on a single line—1.e., the circle |A I. 
A PLANE 
Fic. 2. Course of « and e during variation of point XY VZ — _ = For | 
For the purpose of constructing the solution for a the g 
complete wing, it is found desirable to represent the 
elementary solutions in terms of the complex parameter | 
\ = e”’ and to select a contour C which avoids points on 
the unit circle. It will appear later that the contour 
C can be selected in a way that simplifies the integra 
tion of the elementary solutions over the wing surface. 
In terms of \ the potential of the source becomes | 
g dx 
fe (19) 
47° 1 (ax — By — 2) | 
The integral may now be evaluated by the method of The 
residues. After expressing a, 8, and y in terms of \ alon; 
with the aid of Eqs. (13), we have A SURFACE SHOWING CUT 
ax — By — yz = 
(1/2) [2x — (1 + A*)y — i(1 — A%)z] ~~ (20) TI 
The quantity in the brackets is a quadratic in y and ™ ° 
may be factored so that _ 
: prop 
“-fy- w= - ——= Dr 2-2 (21) | aa 
where 
x—-R yi S58 
6 = —-=* It 
y— 12 x+R pote 
iw t® R_yt+# (22) J 7 then 
“ y—iz x—R theo 
a / x2 "=p Fic. 3. Representation of _— showing two parts of con tegr. 
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This suggests cutting the A surface into two pieces 
aground the unit circle (see Fig. 3). The gap between 
the two portions of the surface can then represent the 
gap between the forward cone of disturbance and the 
rearward cone of disturbance in the (x,y,z) space. If 
the outer portion of the \ surface is now mapped onto a 
second area bounded by a unit circle by means of points 
at inverse radii, the two circular areas can be related to 
the circular cross sections of the Mach cone at x = +1. 
Each point on the interior of the unit circle AI < % 
then corresponds to a point (y + 7iz)/x, defined by a 
ray drawn through the ‘‘Mach circle’ at x = +1. 
The negative ray that pierces the Mach circle at x | 
then corresponds to a point e, on the second surface, 
1/A| < 1. The center of each surface corresponds to 
the center of the Mach cone. When both the fore cone 
and rear cone are considered, it is found that each piece 
of the A surface is covered twice, once by € and again by 
€2. 
Now consider the evaluation of Eq. (19) when the 
contour is drawn just inside the unit circle, enclosing 
1 in the positive direction. For posi- 
+ 3°, the contour will 


the surface |A} < 
tive values of x such that x? > y’ 
enclose €; so that the value of the integral is 
I I | 
27rR 


(25) 


T (y iz) ( € €2) 


For points outside the Mach cone, both «& and € lie in 
the gap between the two portions of the surface and 
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are outside the contour, so that the value of the integral 
is zero. Proceeding toward negative values of x, as 
soon as the point (x,y,z) reaches the upstream Mach 
cone, the pole « moves onto the second portion of the 
\ surface, while the pole €, now appears inside the con- 
tour |A} <1. The value of the integral is now 


al 
27rR 


a (y — 12) (e2 — &) 


(26) 


Exactly the same determination of g, results from a 
contour in the negative direction enclosing the re- 
As later calcula 


mainder of the surface, 1/|y, < 1. 
tions will show, changes in the order of integration and 
differentiation are simplified if the contour C is extended 
around both portions of the A surface and if the points 
\ = Oand 1/\ = 0 (i.e., poles of 86 and y) are excluded 
by small circles. The two parts of the contour C are 
shown in Fig. 3. 

To describe the flow in the region around a lifting 
surface we need an expression for the potential field of a 
‘horseshoe vortex’’ representing the disturbance caused 
by an element of lift at the origin. This expression 
may be obtained by the familiar process of integrating 
the expression for the source in the x direction and then 
These operations 


differentiating in the z direction. 


performed on the integrand of Eq. (19) yield the factor 


—ydd/iax = dB (27) 


so that the expression for the horseshoe vortex becomes 


0, unless x? > y? + 2? 
+xz : 
= — for x > 0 
= (2rR(y? + 27) (28) 


— LS 


2rR(y? + 27) 


for x < 0 


The field is not that of a single horseshoe but of a closely spaced vortex pair extending to infinity in both directions 


along the x axis. 


Outside the cone R = O, the disturbance is zero. 


COMBINED DISTURBANCE FIELD OF A LIFTING SURFACE 


The combined disturbance field for an entire lifting surface s is obtained by superimposing elementary solutions 


of the form (28). 


surface, the strength of the vortices at each point being determined by the local lift /(x,,91). 


This superposition amounts to a double integration of elementary horseshoe vortices over the 


After introducing ap- 


propriate constants and changing the order of integration so that the contour integral is performed last, the ex- 


pression for the combined potential ¢ of the lifting surface becomes 


—] 
Lg ff _ 
Sx’pl ; a(x — x) — Bly — yn) — ¥2 


I(2x1,91) dx; dy, (29) 


It will be shown later that the double integral over the surface s in Eq. (29) yields a two-dimensional complex 


potential function, the form of this function depending on the parameter £. 
By analogy to well-known formulas in two-dimensional potential 


then yields the three-dimensional disturbance. 
theory, the integration over the surface s may be defined 


The final integration over 8 (or) A 


in such a way as to permit differentiation under the in- 


tegral. Performing this differentiation for the velocity components # and @, we obtain from Eq. (29) 








OF 


V S13? pV? 4 5 [a(x 
V Srp V? s [a(x 


The component #/V must, of course, vanish at every 


818 JOURNAL THE 


Relation to Two-Dimensional Flow Theory 


It may now be shown that Eq. (31) represents the 
downwash in the three-dimensional flow by the super- 
position of the downwash of infinitely many two- 
dimensional Each two-dimensional flow is 
associated with an oblique strip drawn in the plane of 
the wing and having its edges tangent to the outline of 
the wing plan form. It will appear that the downwash 
contributed by each strip is given by the familiar 


flows. 


“lifting line’ formula, the loading on each equivalent 
lifting line being obtained by an integration of the sur- 
face loading on the wing in an oblique direction (see Fig. 
4). 

For those parts of the contour consisting of the small 
circles around \ = 0 and 1/X = 0, the equivalent load- 
ing is simply the spanwise loading. It will be evident 
that these parts of the contour yield the vortex drag 
of the lifting surface. Considering first the loop around 
\ = 0, we have 


> 


B (32) 
as — %) — By — 1) — 2 
(1/2d) [ym — (vy + 22)] (33) 


since a(x — x) is negligible by comparison. Eq. (31) 


now becomes 


w@ —] off I(x1,y1) dx; dy, idx (34 
= 34) 
V Srp V? ; ly. — (y + iz)]? A 


LUX,8) = /1ix,¥) a¥ 









i= const. 





WING PLANFORM, S$ 


Integration of surface loading along oblique lines Y 
constant. 
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[(x1,V1) dx, dy, 
: — - a dB (30 
— x) — Bly — m1) — 2]? 
L(x1,¥1) dx, dy, 
y dp (31 


— x) — Bly — ym) — yz)? 


point of the wing surface. 


The integration along x; may now be performed di- 
rectly and results in the spanwise loading, which may be 


f L(x1,91) dx = L’(y1) (35) 


Now the quantity /(x«:,y1) represents the lifting pressure, 


denoted by 


and we may replace it by a quantity that takes equal 
and opposite values on the upper and lower sides of the 
wing the quantity 
L’(y:) may be supposed to take equal and opposite 


surface. Similarly, integrated 
values on the upper and lower sides of the strip, which, 
for \ = 0, coincides with vortex wake of the wing. In 
fact, if we write 

L’'(y1)/pV = AF(y1) (36 
it is clear that AF represents the discontinuity in the 
real part of the two-dimensional complex potential 
function associated with the trailing vortex wake of 
the wing. 


Instead of integrating AF across the strip, we may 
integrate, using the values of F = +(1/2)AF, in the 
positive direction on the upper side and continue in the 
negative direction on the lower side, forming a closed 
contour around the strip. Assuming that the function 
F is smooth and continuous at all points except the 
points of the strip and does not have a pole at infinity, 
the contour around the strip may be deformed so as to 
encircle the point y + iz. Then we have, by Cauchy's 
formula, 


| Jf 1(x1,91) dx, dy, a 
pV Js [vy — (y + iz)? 


F(y1) dy, ; 

: nay = 2miF’(y+iz) (37 

[ym — (y + 22) |? 
The quantity F’(y + iz) is obviously the complex 
velocity function associated with the vortex wake. 
Because of the factor 7, the real part represents the 
downwash and the imaginary part represents the lateral 
velocity v. However, by considering the small loop 
around the point 1/A = 0 on the second portion of the 
\ surface, it is found that the quantity 27iF’(y — 1 
arises instead of Eq. (37). Hence, for the integration 
in the A plane around both loops, the lateral velocity 2, 
which is an odd function of zand discontinuous across the 
strip, vanishes, leaving a real value of the downwash 

w which is continuous throughout the field. 
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Returning to Eq. (34), we have, for the integration 
with respect to the parameter X, 


f idd/d = 20 (38) 


The same value results for each of the small loops 
}—2 0 and 1/A ~ 0. The downwash contributed by 
these portions of the contour then becomes 


w —1 L’(y1) dy, 
= ~ RP. 
| 2rpV* s [nm — (y + 12)]}? 


which is equivalent to Prandtl’s formula for the down- 
wash of the trailing vortex wake. 

It will now be shown that those portions of the con- 
tour C near |\| = | yield the wave drag of the lifting sur- 
face. To illustrate this relation, we replace ax — By — 
yz by the single variable 


(39) 


X = ax — By — 2 (40) 


and for the variable point on the wing we introduce 


X1 = ax; — BY: (4 I ) 
together with the orthogonal variable 
Y, = (Bx + ay) (a? a B") (42) 


The factor a? + £* in the latter expression preserves 


the elementary area. Eq. (31) may now be written 


a 1 off 1(X,,V:) dX, dy, 1g (43) 
= (¢ 3 
Vo Srp? : (X — X;,)? = , 


and the integration with respect to Y; may be performed 
directly by writing 


if wxurn ar. = 0 X.8) 


For \ = | the lines X; = constant correspond to the 
intercepts of a system of plane waves at 45° to the x 
axis (M/ = +/2) and at the angle @ = cos~' 8 in the yz 
plane {see Eqs. (13)]. As the angle @ varies from 0 to 
27, the intersections of the plane waves with the lifting 


= pVAF(X,8) (44) 


surface change their inclination between +45°. For 
these various values of 8, Eq. (44) will correspond to an 
integration of the surface loading of the wing along vari- 
ous oblique directions, as illustrated in Fig. 4. As may 
be seen by introducing Eq. (44) in the surface integral 
in Eq. (43), the downwash contributed by each one of 
these integrated loadings is given by the familiar ‘“‘lift- 
ing line’’ formula 


2ni F’'(X,B) = g — — dX, 
(X — X,)? 


In this form the surface integral in Eq. (31) can be 


(45) 


recognized as the expression for the downwash arising 
irom a given distribution of lift in two-dimensional 


flow. In particular, it is known that, if 


L'(X,,B) = pVAF(X;,B) 


is an ellipse, the downwash, which corresponds to the 
imaginary part of F’(X,8), will be a constant. Hence, 
if the lift is distributed over the wing in such a way that 





> 
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Fic. 5. Elliptic wing. 


the integrated loading in every oblique direction be- 
tween +45° is elliptical, then the downwash contrib- 
uted by the outer parts of the contour C (i.e., |A| > 1) 
will be constant over the entire plan form. If, in addi- 


tion, the spanwise loading is elliptical, then the final 1n- 
tegrated value of the downwash will be a constant. 


DRAG OF ELLIPTIC WINGS AT SUPERSONIC 
SPEED 


MINIMUM 


The foregoing discussion indicated that the downwash 
will be constant over the plan form if the integrated 
This 
condition is, of course, merely a sufficient, and not a 
The condition is met in the case of the 


loading in every oblique direction is elliptical. 


necessary, one. 
elliptic plan form having a uniform surface distribution 
of lift. Hence it is concluded that such a uniform sur- 
face loading yields the minimum drag in the case of the 
elliptic wing. 

The downwash of the elliptic wing may be calculated 
directly from Eq. (31). To evaluate the surface in- 
tegral first, we write 


1(x1,¥1) dx) dy, .. 
, [a(x — x1) — Bly — ym) — vz]? 


2ripVF’'(X,B) (46) 
After introducing / = /) and integrating over x), there 
is obtained 
F'(X,8 l | by | ( lo dyy 7 
(A,B) = —— - 
2rip| ys a\X - ax, + BY: xs 
(47) 


The subscript s has been introduced to denote the limits 
of x1,¥1 corresponding to the edges of the plan form. 
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For the ellipse with semiaxes a and b, we have, in para- 
metric form (see Fig. 5), 


vy = bcos¢; x, = —asing (48) 
and the integral (47) becomes 
lob - sin ¢ d 
F'(X,8) = ——— / —~ -  —_— 
2ripVa Jo X + bBcosd?+aasin ¢d 
(49) 


This form is obviously similar to Eq. (18), the complex 
numbers X, 68, and aa taking the place of the real 
numbers x, y, and z in that equation. 

The evaluation of integrals of this form has been 
given by Jacobi, and a complete discussion will be 
found in reference 10 (see also reference 9). There are 
two determinations depending on the location of the 
poles of the integrand, e; and e2, where 


—X — Vx? — (bB)* — (aa)? 


~~ = 7 
bB — tae 
(50) 
—X + VWxt = (bB)? — (aa)? 
ae = : 
bB — tae 
For real values of 8 = cos 6, the equation 
X? — (bB)? — (aa)? = 0 (51) 


determines those values of X which correspond to planes 
tangent to the edge of the elliptic disc. For points 
(x,y) inside the elliptic disc, 


(x?/a?) + (y?/b?) < 1 (52) 
and in this case we have the general formula 
al cos nd db 
0 X + bBcos¢+ aa sin ¢ 
“oe i sin nd dd 
o X + bBcos¢+ aa sin > 
n n 
a — & ee 
T (53) 


Wx? (bB8)? — (aa)? 














a 
x 
/ 
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id 
4; 
j/ 
/ ; 
ae To 
Fic. 6. Oblique ellipse. 
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Eq. (46) corresponds to n = 


F'(X,8) = lob I : 
Pane eo pVa bB — iaa tof) 


1,* and the formulas give 


For the assumed constant surface loading F’(X,8) js 
thus independent of X, and of x and y, over the surface 
of the ellipse. Hence the downwash will be constant 
at these points. 
tained by introducing the value (54) for (46) in Eq. 


The value of the downwash is ob. 


(31) after making use of the relations 


a=1; y= V1-8° (55 
Thus, 
@ —Iyb V1 — 6? bh of a? 
v sats — mua wie pV? \ .? b? 
(56 
Since the value 8 = 7(a/b) corresponds to two distinct 


poles, one on each portion of the \ surface, the value 
of the integral is twice the residue at this point. 
The drag is now given by 


D = (1 


min V)L (57 
where L is the total lift equal to /,S for a wing of area S, 
The formula for the minimum drag of the elliptic wing 


then becomes 


L? | a* 
Do = : l (58 

, 4(p >) 125 N T pp 
This relation applies at J = +/2. The variation of 
drag with Mach Number may be incorporated in Eq. 
(55) by applying the Prandtl-Glauert rule. The re- 
sulting formula in terms of the coefficients C, and Cp is 


) 1\: 1 \? 
Cres = VM-]1 C,? \()) +( ) (59 
wd 


where A’ = V J? — 1(A) and A is the aspect ratio. 
When the aspect ratio is large, Eq. (59) approaches 


Cp = (VM? — 1/4)C,? (60 


which is the value given by the Ackeret theory for a 
flat wing in two-dimensional flow. On the other hand, 
when A’ is small, the wave drag becomes negligible in 
comparison to the vortex drag, which is given by the 
well-known formula 


C, = C,2/4A (61 


Elliptic Wing at an Angle of Yaw 

Eq. (56) indicates that the drag of an elliptic wing of 
finite aspect ratio is greater than that given by the 
Ackeret theory (for infinite aspect ratio). Smaller 
values of the drag can be obtained, however, by placing 
the wing at an angle of yaw. 


*It will be evident that the formulas for n 
tensions to cases of variable loading over the elliptic wing. 
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MINIMUM DRAG OF 


The treatment of the yawed wing follows the pre- 
ceding analysis with only minor modifications. Again 
the minimum drag occurs when the lift is distributed 
yniformly over the ellipse. With the symbols defined 
as in Fig. 6, the equation of the yawed ellipse becomes 
(62) 


x, = m'y; + (a’/b’)V b”? — yi? 


The poles of the integrand in the equation correspond- 
ing to Eq. (53) now appear, where 

B = am’ + i(a’/d’) (63) 
and there are two distinct poles, one on each portion of 


the A surface. The value of the integral (31) reduces to 


ae ] ( 7 
= R.P. @1 — [m’ 1 (64) 
V pV2S \ + b’ , 
and we obtain 
C,* a’\? 
a - ) i , ‘fel _ 
Crmin 4 mur. y! (m + ¢ «) (65) 


for the minimum drag of the yawed elliptic wing at 
M = v2. 
The variation of drag coefficient with angle of yaw 
vis shown in Fig. 7 for ellipses of various proportions. 
The limit a’/b’ — 0 corresponds to infinite aspect 
ratio, and in this case the expression for the drag co- 
efficient reduces to 


Can. = (C,7/4) RP. V1 — m” (66) 


Minimum Drag Due to Thickness 
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Fic. 7. Minimum drag of elliptic wings at various angles of yaw. 
The change of V1 — m” from a real number to an 


imaginary number as m’ passes through 1 shows the 
disappearance of the drag when the wing 
of infinite aspect ratio is yawed behind the Mach 


wave 


cone. 


To represent the effect of a symmetrical distribution of thickness, we superimpose elementary solutions of the 


form (19), corresponding to sources, over the plan form. 


ties in the combined disturbance field then become 


(1 2)t’ (x1,91) dx, dy, 


u —] g ff ‘ 

V 4r? 5 fa(x — Xi) — Bly - ey yz]? A 
7) | off ‘ 

V Ar? ‘ [a(x — x1) a B(y aun v1) _— yz]? X 


=) 
el 


(1 2)t’(x1,91) dx, dy, 


The expressions for the horizontal and vertical veloci- 


idx 
- (67) 


iy dy 
(08S) 


where ¢’(x1,y1) = Of/Ox, and (1/2)t denotes one-half the thickness of the airfoil. 


In the case of a symmetrical distribution of thick- 
ness, the combined downwash @ vanishes over the wing 
surface, while a value of # remains. The combined 
pressure distribution is given by the relation 


Ap = —puV (69) 


The evaluation of the integral (64) is again especially 
If t’(x1,¥1) 
corresponding to a constant 
Eq. 


simple in the case of the elliptic plan form. 
iS assumed constant 
evaluation of 


The drag therefore 


source density over the surface 
(64) yields a constant value of 7. 
has the minimum value for a given frontal area. Since 





the sections are simple flat-sided wedges, the airfoil 
does not close at the trailing edge, and the figure is a 
semi-infinite body rather than a wing. Since the cal- 
culations are similar to those given for the lifting wing, 


they need not be repeated. The result is 


oe (a?/b*) | (70) 


2 [1/V 1 4+ 


If the source intensity is assumed to vary linearly with 
x, so that ¢”(x1,y:) is a constant, it is found that 0#/Ox 
is constant over the area of the wing. In this case the 
distribution of thickness yields the minimum drag con- 
sistent with a given volume—i.e., 




















JOURNAL OF THE AERONAUTICAL 
(1) 
i LLLLLLL LLL 
12 | 
as 
! L —_ = 
Tr ee a ae 


4 


-8 
La V4 a (2) 
Go / 


, | 
LiL 
2 


0 









































0 / 2 
Axis RATIO b, 


Fic. 8. Minimum drag of elliptic wings having (1) a given 
volume and (2) a given base area. 
" to” | a” -" 
C Dmin. 9 J : r= l sii 2 h2 (4 1) 
a*\vVv 1 + (a?/b?) »” 


In this case the sections have a constant curvature in 
the stream direction. Eqs. (70) and (71) are plotted in 
Fig. 8. 

In each case the area of the cross sections has the 
same distribution along x as that of the corresponding 
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optimum body of revolution. It is interesting to note 
that each of these figures yields the minimum drag for 
all distributions of thickness within the space defined 
by the intersection of its forward and reversed char. 
acteristic envelopes. 
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A Drag Reduction Method for 
Wings of Fixed Plan Form 


E. W. GRAHAM* 
Douglas Aircraft Company, Ine. 


ABSTRACT 


A method is given for reducing the drag of wings of fixed plan 
form, while maintaining constant some property such as wing lift 
or wing volume. This drag reduction is accomplished through 
adjustment of the wing twist and camber if lift is fixed or by ad- 
justment of the thickness distribution if the volume is fixed. As 
an example, the drag due to lift of flat rectangular wings and flat 
delta wings is studied. For the rectangular wings, it was found 
that no drag reduction could be obtained with twist, while cam- 
ber led to a drag reduction that was negligible for high effective 
aspect ratios and large for low effective aspect ratios. The 
camber was most effective near the wing tips. Twist was found 
to be useful in reducing the drag due to lift of delta wings. 
delta wings with supersonic or sonic leading edges, a uniform cam 


For 


ber (constant across the span) produced no improvement. 

The interference drag for two superimposed pressure distribu- 
tions comes from the pressures of one distribution acting on the 
local angles of attack for the other distribution. 
distributions are called orthogonal if their interference drag is 


Two pressure 


zero. The drag due to lift of any pressure distribution may be re- 
duced by combining it with an orthogonal pressure distribution 
This combined distribution can in turn 


A procedure for constructing 


whose lift is not zero 
be improved by a similar process. 


orthogonal distributions is given 


INTRODUCTION 


I’ THE AERODYNAMICS OF WINGS several types of 
problems have been considered. 
mination of the pressure distribution on a wing of given 


is the determination of wing 


One is the deter- 


geometry. A second 
geometry to give a prescribed pressure distribution. 

A third and more difficult problem is the solution for 
the pressure distribution and geometry of the wing 
having minimum drag. This has been investigated 
recently by Jones’? through the use of the concepts of 
reverse flow and ‘‘combined flow fields.’’ He has ob- 
tained criteria for identifying optimum distributions 
and has determined the optimum for elliptical plan 
forms. 

The present paper deals with the third problem using 
an approach different from that of Jones. A procedure 
is found for reducing the drag of a wing of fixed plan 
form by modifying its twist, camber, and thickness dis- 
tribution. By this procedure it is theoretically possible 
to approach arbitrarily near the optimum. 

For simplicity, the basic ideas will be developed in 
connection with a specific case—namely, the reduction 
of drag due to lift for a planar surface of fixed plan 


lorm with supersonic leading edges—though the 
Received July 30, 1952. 
* Aerodynamics Consultant, Santa Monica Division. 


Go 


method can easily be extended to wings having sub- 
sonic leading edges. The problem is thus to find the 
twist, camber, and attitude of the wing which will give 
minimum drag due to lift for a prescribed lift. It 
should be noted that, other conditions being the same, 
the optimum shape of the wing—1.e., its twist and cam- 
will, in general, vary with the amount of lift pre 
However, this variation occurs only in the 


ber 
scribed. 
magnitude of the twist and camber, which will be pro 
portional to the lift. A constructive method will be 
developed so that, starting from any shape of the wing 
which is not optimal, a new shape can be found which 
gives less drag. 

The application of the method to other problems will 


be indicated later. 


DRAG DUvE To LIFT 


We shall consider a wing of fixed plan form at a fixed 
supersonic speed (see Fig. }). All leading edges will be 
assumed to be supersonic. Since only drag due to lift 
will be considered, we may assume, by linearity, that 
the wing has zero thickness. The local angle of attack 
distribution a(x,y), also called downwash distribution, 
will be taken to include both the effect of twist and cam- 
ber and the overall angle of attack. Prescribing a(x,y) 
is then equivalent to prescribing both the shape and the 
attitude of the wing. By the formulas of linearized 
wing theory a(x,y) determines uniquely the pressure 


distribution (load distribution) p(x,y) and conversely a 
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given p(x,y) determines a(x,y). Consider now a dif- Det. = QL1? + doL? (2) mt 
ferent load distribution f(x,y) and its associated angle ! 
Rs ta aga P(x? s § or, replacing L. by Lior. — Li, ott 
of attack distribution a(x,y). If p(x,y) = Cp(x,y), we e the 
shall say that the two distributions are of the same type. Deot. = dyLy? + do(Liot, — L1)? (9 ing 
By linearity, this is equivalent to the relation &(x,y) = _ — ' po 
Ca(x,y). The constant C will be called the intensity of Phe minimum drag for fixed Lto. is found, by differ- 
oe ; : ABs entiation with respect to L;, to be sm 
p(x,y) [with respect to p(x,y)]. Note that for C # 1 I 1, to be tyy 
the two wings have not only different attitudes but also, (Pe) F l ; - 
‘ aes oe ee amet } ae , a = Gy = dop (10 , 
in general, different shapes. Except for flat wings, the a 1 + (d;/de) ia gul 
wing corresponding to &(x,y) is not obtained by a rigid ; : aim 
rotation of the wing corresponding to a(x,y). If p(x,y) and is obtained when for 
is not a multiple of P(x), the two wings (or their load L1/Ltot. = d2/(di + de) (11 Th 
distributions) will be said to be of different types or in- my — as} 
dependent Since d; and d, are always positive, it follows that . 
Consider now two independent pressure distributions dopt. is less than d, and less than dy. So, having given pr 
pi(x,y) and p»(x,y) with corresponding downwash dis- @" arbitrary pressure distribution, if an orthogonal dis. a 
tributions a;(x,y) and a:(x,y). If a wing has the first een (with finite d) a then the drag at on 
type of pressure distribution with intensity C,, then its fixed lift can be reduced by shifting the correct amount - 
rpm of lift into the orthogonal distribution. This is true | 
even when d for the orthogonal loading is worse than 
' , : — ‘ ae edg 
D, = f Cipi(x,y) Craa(x,y) ds (1) d for the original loading. However, if the orthogonal a 
loading carries zero lift, it will have an infinite d and ‘ f 
i e ° ° ° ° ° tae ’ 2 e er 
If the two types of distribution with intensities C; cannot be used for drag reduction. - 
: ! for 
‘ > are ¢ > sec > al drag 1s _— ° ° ° een 
ind C; are superimposed, the total drag 1 rhe preceding development indicates the possibility is | 
 —— f [Cipilx,y) + Copo(x,y)] X of drag reduction after certain orthogonal loadings tro 
: s " , : : have been found. The procedure for finding such not 
[Cyai(x,v) + Crae(x,y)] ds (2) a aes a ; big se Te 
- : loadings is as follows: Given an arbitrary distribution val 
Regrouping the terms Pil(x,y), ai(x,y). Let po(x,y) be defined by I 


P re 12 Copo(x,y) = C,py(x,y Py (Xr 2) ia 
Dot, — f [C 1P1a%1 =f c 3Pr2a | ds + 2P2(Xx - ) CuPa(x y) + CyPo | ¥,V) (J = 
CC. f [Pras + frm] ds (3) where p, and p, are arbitrary but independent types of 


loading. Then, if p2 is to be orthogonal to /, the fol- 


The second integral corresponds to the interference lowing equation must be satisfied: 
drag of the two loadings. A pair of ‘orthogonal’ a : ; : 4 . ; , 
loadings will be defined by setting this interference f [(Caba + Cobr)Ciar + (Caaq + Cyay)Cipi] ds = 0 = 

9 ‘as 

drag = 0. (13 ¥ 
con 

Thic tte ‘ 2c > Narace- y are > 7 ? bd P 1, 
f Ibis + Prd ds = 0 (4) his determines the necessary value of C,/C,, and an ber 
s orthogonal distribution p, is given by of : 
Orthogonality is then a property of the ¢ypes of pressure ° F for 
distribution and is independent of their intensities. J. (Pia, + pao) ds dra 
’ si ite aia ” p2 = Py — Pr (14 
_Assume that p; and p, satisly t aS orthogona uty con- f (pia, + par) ds I 
dition. Then the drag of the loading Ci/i(x,y) by itself P clu 
or in combination with Cypo(x,y) is If p, is chosen to be identical to /p,, then son 
I 
Dd, = cef pia ds (5) 2. pres ds . the 
p: = pi - ; Pr (15 | 
and its lift is f (pia, + pyar) ds : 
§ or 

In = G J pi ds (6) Thus if fp; is a given pressure distribution and #, is by 
an arbitrary but independent distribution, then an rele 
the ratio D,/L,’ is given by orthogonal loading pf. can be obtained. Unless this | 
° 2 = loading carries zero lift, it can be used to obtain a drag rec 

dD, /L;? => f Pia ds (f Pi ds) (4) 6 _ 

s s reduction. 

Since this ratio is independent of C, (the intensity of This type of procedure can be repeated to obtain a Mo 
loading), it is a convenient characteristic of the type loading ~; which is orthogonal to p; and p.. The con- 
of loading. To simplify the notation let d, = D,/L,°.  tinuation of this process makes it possible to approach 


For the pair of orthogonal loadings p; and ps, the arbitrarily close to the optimum, which will contain all 
total drag is the lifting members of the orthogonal set. The optt- 
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DRAG REDUCTION FOR 
mum must be orthogonal only to zero-lift loadings, 
otherwise it could be improved. It can also be shown 
that the optimum is orthogonal to every zero-lift load- 
ing. 

In practice it may be sufficient to consider only a 
small group of orthogonal loadings. Fig. 2 shows 
typical drag reductions obtained by using a single load- 
ing orthogonal to the flat-plate distribution. Rectan- 
gular and triangular plan forms are used for these ex- 
amples, and in some cases alternative a, distributions 
for determining the orthogonal loading are included. 
The greatest improvements are shown for low effective 
aspect ratios. 

For the rectangular plan form, twist offered no im- 
provement, but symmetrical camber can be used to re- 
duce the drag. This is also true for any plan form sym- 
metrical about an axis (in the plane of the plan form) 
perpendicular to the flow direction.’ 

For the delta wing with sonic or supersonic leading 
edges, a simple camber produced by bending the wing 
along spanwise lines gives no reduction in drag. (In 
reference 3 it is shown that a similar conclusion holds 
for a more general class of plan forms.) However, it 
is possible to improve the flat-plate delta wing by in- 
troducing twist. It is probable that a camber that is 
not uniform along the span could also be used to ad- 
vantage. 

For more general plan forms, both twist and camber 
would be required to attain the optimum shape for 


lowest drag due to lift. 


FURTHER DEVELOPMENTS 


The method used for reducing drag due to lift can 
also be applied to reducing the thickness drag. In such 
cases the wing volume or the frontal area might be held 
constant, or the thickness required for structural mem- 
bers could be specified. The significant characteristic 
of a type of pressure distribution would then become, 
for example, drag/(volume)* in the first case, instead of 
drag/(lift)?. 

If the definition of orthogonality is extended to in- 
clude leading-edge suction, then plan forms with sub- 
sonic leading edges can be studied. 

It is also possible to use similar methods in reducing 
the drag of wings in the upwash field of a fuselage. 

A number of mathematical theorems relating to 
orthogonality and reverse flow have been developed 
by A. M. Rodriguez* and will appear in the report of 
reference 3. 

A study of the drag reduction problem for delta and 
rectangular plan forms is currently being made by Miss 





* Aerodynamicist, Douglas Aircraft Company, Inc., Santa 


Monica Division. 
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Drag reduction at constant lift. 
Beverly Beanet and Kelsey Walker.t The examples 
they have considered will appear in the reports of refer- 
ences 4 and 5. 

It should be mentioned that drag reduction methods 
can be constructed without reference to the concept of 
orthogonality. For example, this has been done by 
Adams and Sears. The method used in the present 
paper was chosen because it seemed best adapted to a 
general analysis of the problem. 
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On Curved Shocks in Steady Plane Flow 
of an Ideal Fluid’ 


C. TRUESDELLt 


Indiana University 


SUMMARY 


For an inviscid fluid obeying an arbitrary equation of state, the 
jump in the velocity gradient across a curved shock in steady 
The results are expressed in terms of 


plane flow is calculated. 
For the case of uniform oncoming flow, 


dimensionless variables. 
it is shown that the vorticity generated by the shock is a function 
of curvature, strength, and shock angle only, being independent of 
the form of the equation of state. 


INTRODUCTION 


| eameangl SHOWED THAT the discontinuity in the 
velocity gradient at a shock wave is expressible in 
terms of the discontinuity in the velocity. Thomas? 
has introduced certain matrices that make practicable 
the explicit calculations for the case of steady plane 
flow of a perfect gas with constant specific heats, de- 
void of viscosity and thermal conduction and subject 
to no extraneous force. Our purpose here is to extend 
these results to the case of a fluid obeying an arbitrary 
equation of state, at the same time putting them into a 
more compact form in terms of dimensionless variables. 
Since the method is precisely that of Thomas, we do not 


explain the details. 


GENERAL CASE 


Quantities written with a subscript 1 are evaluated 
at the front side of the shock; quantities without sub- 
scripts, behind it. Then, in the usual notations, the 
Rankine-Hugoniot conditions assume the form 


Pilly, = Puy, 


Pi + Pilly,” _ p = pu,” 


(1) 
Uy = U; 
hy + (1/2)u,,.? = he + (1/2)u,,? 
We introduce the following abbreviations: 
[p| u u 
6 => P 7, => -, x >= : 
Pl Uin Uy 
(2) 
u u, 
Mi», S , M,, = 
Cy Cc 


where 6 is the strength and 7 the obliquity of the shock. 
Denoting the rectangular coordinates of velocity by u, 
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* This work was done under Navy Contract No. N6onr-180, 
Task Order V, with Indiana University. 

+ Professor, Graduate Institute for Applied Mathematics. 


826 


and of the unit normal to the shock line, directed toward 
the region behind the shock, by v,, we have then, from 
Eq. (1), 


+s [ue] OU Va 
p = ( 5p, [u.] = — 
Pl re 

Un 1 +. 5 
Uy, = ’ = )T (3) 

ita * ( 

611 + (1/2)6} 

bl = U,,°, [h| = Sl,," 
(P) = jr Ih = a 
Since h = h(p, p), a formula for 6 as a function of u,,, 


pi, pi Or M,, pi, pr: may be obtained in principle by 
elimination among Eqs. (3)s, (3)5, and (3);. For the 
special case of a perfect gas with constant specific heats, 


we have 


8 = (M,,? — 1)/[1 + (1/2) (y — 1I)M,,7] 
Urn? l 6 
M.° = pi : = + 
Y~P\ 1 — (1/2) (y — 16 
1— M,? = (4) 
a + 1 M,,” or l 
- ) —_ 
3° yM,,? — (1/2) (vy — 1) 


(1/2) (vy + 1)6 
l + (1/2) Wy + )6 


In addition to the usual equations of motion, 


Palla + Pug a _ 0 / = 
0 | 


Pa = Puguag 
we require the following form of the condition that the 


entropy is constant along each stream line: 


Ue ggg — CU, = 0 (6) 


a,f 


To derive it, note that, since n = n(p, p), we have 


= Ue Ue _ ) alla Palla 
Op p Op p 


whence, by Eqs. (5), follows 


(On/Op) » 
(On/OP), 


Ug glial g + co = O 
but c? = (0p/0p),,. 

Let a prime denote differentiation with respect to 
arc-length defined so as to increase algebraically in the 
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CURVED SHOCKS IN 


direction of the unit tangent vector A, = @gg ¥g- From 
Eqs. (3)2, (3)s, and (2), follow Thomas’ equations: 

Ug = Ua prtg = tyept\g + Ag=A.* 

pb’ = Padrtg = Pipts + B= B* (7) 
Pi,a\g +C=C* 


( ‘) 
Uin” | » 
i-a" 


C = (p,6)’ (8) 


/ 


p _ P.ad~ ai 


( stint) 
“ ; B 
1+ 6 


Note that the obliquity 7 does not enter these formulas 
explicitly. The complexity of the results can be re- 
duced by the abbreviation 


where 


l 


‘ll 


= (1 + 8) M,,(06/0M,,) (9) 


but we shall not go into the details. For a perfect gas 
with constant specific heats, Eq. (4) put into Eq. (9) 


yields 
4 — | 
4m (1 a i) (10) 
a ee 2 ; 
Now set 
Cu Ca Ar As 
= | Uu uo 
Cr Cx 

Uy, Uy, 
Dy Da U2 Vi} 
Dy» Dx = = = v2 

Uy» 


Then DjaCa; = 6;;. Note also that 


1 t = 
Dll =| 7 * ~% (11) 


D 
see 


la 


Let quantities B;,; be defined as follows: 
Ku,B;; => Ue, pms Bj (12) 


where x is the curvature of the shock line. This 
definition fails when x = 0, a possibility we exclude 


henceforth. From Eqs. (7) and (6) we obtain 


B* 
si, - = 
ss By By _ pu, (13) 
“"" |\|Ba Bell ||\Ao*te ie « 7 
tn M,,’ 
Inversion of Eq. (12) yields 
Usp = KUyB,;Dj.Djg (14) 


From Eq. (11), then, 


Ugg = ku, [Bu(l + x?) — 2Baynx + Be] (15) 
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whence, by Eq. (13), 


a By(l + x?) — 2Bazx 


(16) 
” M,? —1 ; 


Thus B,: like the other B;, may be calculated without 
knowledge of u,3, so that Eq. (14) gives the solution 
of our problem. 

To complete the analysis, note that by Eqs. (7) and 
(5): we have p,Cs; = a; where a = C*, w = 
—pkM,,? Bw. Hence, 


Pa = Dia (17) 
From Eq. (5): follows 
Pa — pu,” BiDia (18) 


For the squared magnitude of the pressure gradient we 
have, by Eq. (11), 


P a? a = PKU, BieBjeDiaD ja 


= p*x*t.,* [ By»? 1 x By — By)? | (18A) 


and, hence, the theorem that the pressure gradient be- 


hind a curved shock vanishes if and only if By = 0, 
By = 0. 
The geometrical formulas \,’ = kvg, Yq’ = —KAq 


enable us to calculate the derivatives of various quan- 
tities along the shock line. 


u ‘= KU,Bi Dia 


a 


Un! = k(u,,.Boy — u,(1 + By) | 

(19) 
u,’ = ku,By, ve’ = xu,7Bo; 
p’ = oy p’ _ —kpu,,"By» 


while Thomas’ formula for the curvature A of the 
stream lines behind the shock becomes 


—_ Bus 
K«—— (—B: + ~= ) (20) 
V1l+ x? Se & 


Eqs. (19)s, (19)4, (19)6, and (13), give physical mean- 
ings for the quantities B;;. 


Let w be the vorticity. Then, 
w= 2xu, Bro (21) 
Hence, 
wu, = vw’ + (p’/p) (22) 


CASE OF UNIFORM FLOW AHEAD OF THE SHOCK 


When pressure, density, and velocity are uniform in 
the region in front of the shock, not only is A,* = 
A,, B* = B, C* = C, but also the differentiations (8) 
are easy to carry out with the aid of the relation u,,’ = 
expressions for the B,; are 


—xu, The resulting 








Bu 
Bn Bell ~ ||By + 78? 


! ' 





In any fluid, shock strength 6 is to be expected to be an 
increasing function of M,, in all circumstances. It will 
then follow by Eq. (9) that A 2 0. From Eq. (10) it is 
plain, since 6 S 2/(y — 1), that such indeed is the case 
for perfect gases. By Eq. (23) we conclude then By, > 
0, Ba > 0, provided only r + 0. Comparing this 
fact with the italicized statement following Eq. (18A), 
we obtain the theorem that when a uniform flow of any 
fluid breaks across a shock, the pressure gradient cannot 
vanish on the rear side of the shock at any point where 
the shock is curved and oblique. Also, 


C = —pirxA(1 + 3) (24) 





For a perfect gas with constant specific heats (23) can be 
| 38 


| .|| = 
Bi | By + 76? 


while, for Eq. (24), we have 


— 
7? 3 





ates ie (1-757 s)a+o (28) 


9 


~ 


Substitution of these special values into the general 


[y +1 — 7°11 + 1/6) (8 + [7 — y]8)] [1 + (1/2)(y + 16] 


r(A + 26) 
6 — 7r2(1 + 6) (2A + 36) (23) 
M,’? - 1 
For the vorticity we obtain the remarkable formula 
Uy ,7 8" u,b" 
w= Ku,7Td? = a x (25) 
1+ 6 1+ 6 


showing that the magnitude of the vorticity generated by 
a shock of given strength and curvature depends only on 
the magnitude of the tangential component of velocity and 
is independent of the form of the equation of state. An- 
other remarkable property of the vorticity follows by 
comparing Eqs. (23), (19)3, and (25): 


w= u,'7rd (26) 


expressed explicitly in terms of 7 and 6 only. 


4r(1 + 6)/(1 + vy) | 


- (27) 
(1/2) (y + 1) 


formulas, Eqs. (14) and (17)-(20), yields the various 
results of Thomas. 
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Elastic Stability of Simply Supported Flat 
Rectangular Plates Under Critical Combina- 
tions of Longitudinal Bending, Longitudinal 

Compression, and Lateral Compression 


ROBERT G. NOEL* 


North American Aviation, Inc. 


ABSTRACT 


A theoretical investigation was made of the buckling of simply 
supported flat rectangular plates under critical combinations of 
bending, longitudinal compression, and _ lateral 
Interaction curves for these loadings are presented 


Also included are curves of 


longitudinal 
compression. 
for various plate aspect ratios. 
critical buckling constants for simply supported flat rectangular 
plates under combined unsymmetrical bending and lateral com- 
pression. The results indicate that the reduction in the allow- 
able bending stress due to the addition of lateral compression is 
greatly magnified by the further addition of only a small longi- 
tudinal compressive load. 


INTRODUCTION 


y THE DESIGN OF STRESSED-SKIN STRUCTURES, con- 
sideration must sometimes be given to the critical 
stress of a plate under combined longitudinal bending, 
longitudinal compression, lateral compression. 
Box-beams webs, for example, are subjected to these 
loading conditions due to the bending moments alone. 
Also in beams having unsymmetrical distributions of 
cross-sectional area, the webs are subjected to unsym- 
metrical longitudinal bending stresses and lateral com- 
pression (crushing) stresses due to the applied moment. 
In the present paper, the critical stress coefficients 
are found for the plate loading shown in Fig. la. 
Critical stress coefficients have been calculated by 
other investigators for certain boundary values of the 
more general loading case assumed here. Some of these 
are: longitudinal bending and compression,' Fig. 1b; 
longitudinal bending and lateral compression,’ Fig. 1c; 
and compression in two directions,*» * Fig. ld. The 
results of these papers were used to aid in the con- 
struction of the interaction curves presented herein. 


and 


SYMBOLS 


a = Jength of plate in x-direction 

a, = Fourier coefficients 

b = width of plate in y-direction 

D_ = plate flexural rigidity, Et?/12(1 — yu?) 

E = Young’s modulus of elasticity 

k, = critical longitudinal bending stress coefficient 

kez = critical stress coefficient for longitudinal bending and 


compression combined 


Received May 26, 1952. 
*Stress Analyst, Aerophysics Laboratory. 


829 























































kz = critical longitudinal compressive stress coefficient 
k, = critical lateral compressive stress coefficient 
m = number of half waves in x-direction 
N, = longitudinal bending force per unit length acting on 
plate middle plane 
No = maximum compressive value of VN; at edge of plate 
N, = longidutinal compressive force per unit length acting on 
plate middle plane 
N, = lateral compressive force per unit length acting on plate 
middle plane 
No: = sum of longitudinal loadings N, and Ny 
R_ = stress ratio 
T = external work of applied forces 
t = thickness of plate 
V = internal energy of deformation 
w = deflection of buckled plate normal to middle plane 
y 
Nbx * No(I - oc —) 
x 
- y 
Nbx } GENERAL 
LOADING 
y Ny 
(a) TYPE OF LOADING SOLVED IN PRESENT PAPER 
a “ad 
meee Ny =O 
\— = “N y = 
g bx 
a —— 
(b) TYPE OF LOADING SOLVED IN REFERENCE |. 
Ny 
bbb bbb db bbsd dtd 
= bh ; 
y 7 o*2 
a = 
—— =" 
Tht} tI ttt) 
(C) TYPE OF LOADING SOLVED IN REFERENCE 2 
Ny 
—peeee bb bbe bb 
a ~ 
— = Nx _— 
TTT TTT 
(d) TYPE OF LOADING SOLVED IN REFERENCE 3 
FIG.1 BUCKLING OF A PLATE UNDER COMBINED LOADS. 
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fe) 04 08 1.2 16 20 24 28 32 
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FIG 2 CRITICAL STRESS COEFFICIENTS FOR FLAT 
RECTANGULAR PLATES UNDER COMBINED LOADS 


x,y = plate direction coordinates 

y = distance from compression side of plate to the neutral 
axis in bending 

a = ratio defining location of axis of zero stress b/y 

ratio of wave length A to plate width 6, = »/b = a/mb 


~ @® 
ll 


ll 


half-wave length of plate buckles in x-direction 
= Poisson’s ratio 


= stress = N/t 


THEORETICAL ANALYSIS 


The critical stress coefficients for an infinitely long 
flat under the loading of Fig. la is found by the well- 
known energy method (Rayleigh-Ritz). 

Consider the simply supported plate loaded in the 
middle plane by the evenly distributed lateral forces 
N, and the longitudinal distribution given 
by 


force 


Noz = Noll — a(y/b)] (1) 


Then various longitudinal force distributions can be 
obtained by changing a. If the plate is assumed to 
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buckle into a sine wave in the x-direction, the deflected 
surface w can be taken in the form 


w = sin (rx/)) >. a, sin (nmy/b) (2) 


Because each term of Eq. (2) and its second derivatives 
vanish at x = 0, x = a and at y = 0, y = 3, it is clear 
that the boundary conditions of zero deflection and 
moment are satisfied at the edges of the simply sup- 
ported plate. 


The elastic-strain energy of bending is given by 
Eq. (199) of reference 1, 


y= D f° B= 4 a) 21 ." 
2 Jo 0 Ox? Oy? 7 ? 
Ow Ow o’w \* || 
eer -dxdy (3) 
Ox? Oy? dx dy/ J 
By substituting Eq. (2) and integrating, 


Dybx* * =" l n*\? 
V = : Pig t) 
Ss >» (<. 7 *) (4 


The work done during buckling by the forces acting in 
the middle plane of the plate is given by Eq. (201) of 


‘fas dy 


reference 1, 


] b r ‘0°w\" O*w 
ToL Ll Ga) +S) 
2 0 0 Ox? oy? 


(9) 
By substituting Eqs. (1) and (2) and integrating, 
Nobr? ( “) _=* 
T = 1 —  an*+ 
16\ Y 
Nae” = asim gWMe eee ge 
>: > = ao : > n-a," (06) 
d n=l ¢ (9° — 8°)° Sb n=1 


(n + i must be odd). By equating the work 7 to the 
internal strain energy V during buckling and making 


the substitutions 
No = k,,3*D b? 


N, = k,w?D/b? 
B /b 


an expression in k,, and k, is obtained which is to be 
minimized by adjusting the a,, coefficients. 


n U 


> ee 
> a,” [1 — n°6?]? =k,Bt >> n2a,2+ 
n=1 n=1 
Pa igo Sk,,B?a "=< a,amnt 
kyr B? (: - *) > a,?+ — 2 2 
“ n=l Tv n=1 a a ET 
(7 


(n + 4 must be odd). For a given k, constant, the ex- 
pression for k,, is minimized with respect to a, by 
equating to zero the first derivative of Eq. (7) with 
respect to a,. On simplifying, this gives the following 


set of linear homogeneous equations: 
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Qa 
t E — n°B)? — kyBin® — ky sB (: 3) 7 


Sk,:B?a = 0 
7 a (n* — $*)* 


ant 
i (s)* 


(n + 7 must be odd). For purposes of calculations, it 
was found that sufficient accuracy is obtained by 
setting : 


Qj+6 = 0 ¢ _ 0, 2, 4, cag: SS 


In this problem, it is also convenient to solve first for 
the odd coefficients in terms Of the even coefficients 


Sk,,B?a 2a2 da, 6d¢ 
: = kore Se (9) 
wA (727 — 4)? (22 — 16)? (2? — 36)? 


@=1,3,5,7¢,... 


aj = 
o ), where 
A = {(1 + #8)? — k,B4? — ky.6*[1 — (a/2)]} 


On substituting Eq. (9) into Eq. (8) and simplifying, 
the following set of stability equations is obtained: 


a 7 
On E + np")? — k,n?B4 — k,,B? (: > *) - 


64k,,"Bta°n > 72 
ania Sea = ae — 
m4 7 (n? — 1?)2(42 — 4)2A 
64k, ,*Biarn 2 172 
> SSS (10) 
1 ; (n* — i7)?(42 — 16)?A 


where (1 = 0, 2,4: 4 = 1,3, 5,..., @) and 


A = {(1 + 7°67) — k,i*84 — kyB?[1 — (a/2)]} 


The only nontrivial solution of Eq. (10) is found by 
setting the determinant of the coefficients equal to zero. 
For any given k, constant, the lowest value of k», is 
the desired value of the buckling-stress coefficient. 
Having first evaluated the infinite summations, it was 
only necessary to find a sufficient range of values of ky, 
and k, to satisfy the above conditions. 


RESULTS AND DISCUSSION 


a) Effect of Finite Plate Length 

To obtain the results of Figs. 2, it was only necessary 
to make the substitution B = a/mb and choose the 
value of m such that k,, is a minimum. The values 
obtained by this substitution are valid for all plate 
a/b ratios where the length a is divided into m half 
waves of length X. 


(b) Extension of Result Curves 


In Figs. 2, the results of the above calculations are 
given for values of a = 7/4, */s,5/s, and 1. Sufficient 


calculations have been made so that values of ,, may 


* Note that this set of equations reduce to Eqs. (7) reference 2 
by setting a = 2 and are equal to Eqs. (g), pp. 353 of reference 1 


for k, = 0. 
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FIG. 2 CONTINUED 


be found for any k,, 8, and @ by cross plotting the re- 


sults of the present paper and references 2, 3, 
and 4. 
It can be seen from Eq. (1) that, for a value of a = 2, 


the case of pure bending is obtained,’ and setting a = 0 
gives the case of uniformly distributed compression in 
two directions.*4 


(c) Interaction Ratios 


The use of stress ratios sometimes provides a con- 
venient method of representing critical stresses of plates 
under combined loads. The ratios of the stress actually 
present in the plate to the critical value when no other 
loads are present is called the stress ratio. 

Using Figs. 2 of this paper and the results of refer- 
ences 2, 3, and 4, the interaction stress-ratio curves of 
Figs. 3 were constructed. Because the stress ratios 
depend to a large extent on the plate aspect ratio, the 
curves are given for values of a/b = 0.8, 1.0, 1.2, 1.6, 
2.0, 3.0, and ~. 
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INTERACTION STRESS RATIOS FOR SIMPLY 
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SUPPORTED 


FLAT PLATES UNDER COMBINED LOADING. 





Critical stress coefficients for simply supported plates 
under each loading N,, N,, and N, are given in Fig. 4 
for reference. 

Also, it was found that the combined loadings of 
this report result in critical stresses that cannot be de- 
scribed by the usual nondimensional equation of the 
type, 


R,? + R.! + Ri + on™ l 


because of the dependency on the plate aspect ratio. 


USE OF THE CURVES 


(1) The result curves can be used to design a plate 
to resist buckling under given loads N,, and N,. The 
buckling constants are first established for any plate 
aspect ratio from Figs. 3 by selecting a point for which 
the ratio of the k values is equal to the ratio of the 


applied loads. 


k, Ro, = Ne No 
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FIG. 3 CONTINUED 


1.0 


08 


06 


04 


O02 











By using the k values that satisfy this condition, the 
critical skin thickness can be found by solving either 
of the following equations: 

N, = kywEt*®/12(1 — w?)b? 

No = Ror Et®/12(1 — w?)b? 
By making the thickness ¢ slightly greater than the 
computed critical value, a positive margin of safety is 


obtained. 


(2) In order to determine if a plate of a given 
thickness is stable or unstable under given loads, the 
use of the stress ratio curves of Figs. 3 is more con- 
venient. First read the ky values from Fig. 4 and calcu- 
late the critical stresses for each loading alone. 


Tx = k,, [rEt?/12(1 — w*)b?] 
Ty ind k,, [r?Et? 12(1 —_= u”)b?| 


o», = ky, [w?Et?/12(1 — w*)d?] 
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FIG. 4 CRITICAL STRESS COEFFICIENTS FOR SIMPLY SUPPORTED FLAT 


PLATES UNDER A SINGLE 


Next calculate the applied stress ratios R’. 
R,’ = (N,/t)/o2 
R,’ = (N,/t)/ey 
R,’ = (N, t) Tv 


II 


Plot two of the applied ratios R,’ and R,’ on Fig. 3 
The third applied ratio R,_ given by the curves at this 
The stability of the 
given plate under the applied loads is defined by com- 
paring the critical load R 


point is the critical value of R,. 


r,, With the applied value R,’ 


R,,>R,' (stable) 
R,. = =i (marginal) 
R,,< R, (unstable) 


LOADING —-— Ner=k 1*D/b* 
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Lateral Instability of Rectangular Beams of 
Strain-Hardening Material Under 
Uniform Bending 


W. H. WITTRICK* 
Unwersity of Sydney, Australia 


SUMMARY 


The problem considered in this paper is that of the lateral buck- 
ling of initially straight beams of thin deep rectangular cross sec- 
tion under uniform bending, when the bending moment is large 
enough to cause plastic flow to occur over part of the cross section. 
The solution applies to any shape of stress-strain curve and repre- 
sents a generalization of that given by Neal! for the case of mild 
steel bars. 

Numerical calculations are made for two different materials 
one representative of an aluminium alloy and the other of an- 
nealed mild steel—and the two materials are compared. 


(1) INTRODUCTION 


I A THIN DEEP BEAM that is initially perfectly straight 
is subjected to a bending moment about its major 
principal axis, collapse may occur by lateral deflection 
and twist occurring simultaneously. Solutions to the 
problem of determining the critical bending moment 
when the beam remains elastic have existed for many 
years, and recently Neal! has carried out theoretical 
and experimental work on the same problem for mild 
steel beams of rectangular cross section in which part 
of the section yields before lateral buckling occurs. 
The present paper extends this work to materials in 
which strain hardening occurs and is applicable to any 
shape of stress-strain curve. 

In the classical theory of elastic stability it is as- 
sumed that, when the critical load is reached, buckling 
occurs while the load remains constant. It is known, 
however, that in problems of instability in the plastic 
range there is a somewhat lower load above which lat- 
eral deflections can and will occur with increasing load. 
In the present paper these two loads have been called 
the upper and lower critical loads, and they mark the 
boundaries of a region in which buckling takes place. 
Being guided by experimental results for the plastic 
buckling of axially compressed struts, we assume that 
the collapse load lies somewhere in this region. 

Neal’s experimental work on lateral buckling of rec- 
tangular beams under uniform primary bending mo- 
ment was confined to the case of annealed mild steel 
that possesses the phenomenon of an upper and a lower 
yield stress. His results showed extremely good agree- 
ment with the upper critical bending moments as cal- 
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culated theoretically. However, it is shown in this paper 
that, for the particular material that Neal used, in 
which the ratio of the upper to the lower yield stress 
was 1.4, the difference between the upper and lower 
critical bending moments is small, of the order of 0.5 
per cent. It is therefore felt that these experiments 
do not provide conclusive evidence that the collapse 
load is given more accurately by the upper critical 
bending moment than the lower one. 

Calculations for a material having a stress-strain 
curve representative of an aluminium alloy are also 
described in this paper. In this case it is found that 
the difference between the upper and lower critical 
bending moments is of the order of 5 per cent, which is 
considerably greater than for the annealed mild steel. 


The analysis shows that it is far easier to calculate 
the lower critical bending moment than the upper one, 
and it is felt that this will also be true for 
other, more complicated, sections than for the 
simple rectangular one considered here. Also, the 
lower critical bending moment should give a slightly 
conservative estimate of the collapse load, and 
for these reasons it is suggested that it is probably of 
greater practical importance than the upper critical 
bending moment. 


(2) CONDITION FOR LATERAL DEFLECTIONS TO OCCUR 


The problem considered in this paper is the lateral 
instability of a perfectly straight beam of narrow rec- 
tangular cross section when subjected to terminally 
applied bending moments about its stronger, or major, 
principal axis. If the material remains elastic, when 
this primary bending moment reaches a certain critical 
value the beam deflects laterally and at the same time 
twists about its longitudinal axis, the two end sections 
being prevented from rotating. Collapse then occurs 
under constant load. The expression usually given in 
textbooks (see, for example, Timoshenko’) for the 
critical bending moment /,, can be written in the form 


1/L = (M.,/®) V1/BC (1) 
where L is the length of the beam; B is the secondary 


flexural rigidity for bending about the weaker, or minor, 
principal axis; and C is the torsional rigidity. 
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In the derivation of Eq. (1), the curvature of the 
beam due to the primary bending moment is neglected 
entirely. This is justifiable only if A, the primary flex- 
ural rigidity for bending about the major principal 
axis, is extremely large compared with B and C—i.e., 
when the ratio of the width to the depth of the cross 
section is small. A more accurate discussion of the 
problem was given by Reissner® for the case of a canti- 
lever beam subjected to a shear force applied at the 
centroid of the free end in a direction producing bending 
about the major principal axis. He used Kirchhoff’s 
equations for the combined bending and twisting of a 
rod (see Love‘) to take account of the curvature due to 
the primary bending moment, and Neal! has adapted 
this method to obtain the critical bending moment in 
the case of terminally applied couples. The formula 
that Neal obtains is 


1 Me ( x “C a ') (2) 

Lo YW\B A/\C A " 
and it will be seen that this reduces to Eq. (1) if the 
primary flexural rigidity A is extremely large. 

If plastic flow has occurred under the action of the 
primary bending moment before lateral deflection takes 
place, it is necessary to decide exactly what is meant by 
the critical bending moment. In principle, the prob- 
lem is similar to that of the buckling of a straight strut 
under axial compression when the buckling stress ex- 
ceeds the elastic limit of the material. If it is assumed, 
as in the corresponding elastic problem, that the strut 
remains straight until a certain critical load is reached 
and that buckling then occurs under constant load, the 
value of this load is given by the reduced-modulus or 
double-modulus formula of von Karman. It is shown 
by Shanley,’ however, that this is not the smallest load 
at which the strut can assume a bent equilibrium con- 
figuration. In fact, lateral deflections can occur with 
increasing load at a value given by the tangent-modulus 
formula of Engesser. Shanley also concludes that the 
collapse load—i.e., the maximum load which the strut 
can carry—lies somewhere between the tangent-modu- 


lus and reduced-modulus loads. Experimental results 
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seem to indicate that the collapse load is usually much 
closer to the tangent-modulus than the reduced- 
modulus load, except if the strut is short, in which case 
the reverse is often true. Thus, in the plastic range, it 
is not possible to give a clear-cut definition to the criti- 
cal load. We can, however, define the lower critica] 
load as the load at which lateral deflection can occur 
with increasing load and the upper critical load as the 
load at which lateral deflection can occur under con- 
stant load. 

In the problem discussed in this paper, it is evident, 
by analogy, that there are two similar extreme cases to 
consider: (a) to determine the lower critical bending 
moment at which lateral deflection and twist can occur 
with increasing bending moment; and (b) to determine 
the upper critical bending moment at which, if no pre- 
vious lateral deflection had occurred, the beam would 
buckle laterally under constant bending moment. 

By a further analogy with the strut problem, it will 
be assumed that the maximum bending moment that the 
beam can carry lies somewhere between these upper 
and lower critical values. 

Neal! shows that, if more precise definitions are 
given to the rigidities A, B, and C, Eq. (2) still applies 
in the plastic range. Thus, A is now defined as the 
overall primary flexural rigidity—1.e., the ratio be- 
tween the primary bending moment and primary curva- 
ture immediately before lateral deflection occurs. Its 
value obviously decreases as the amount of plastic 
flow due to the primary bending moment increases. 

The term B, which Neal calls the initial secondary 
flexural rigidity, is defined as the ratio between an in- 
crement of secondary bending moment, applied about 
the minor principal axis, and the corresponding incre- 
ment of secondary curvature. For the determination 
of the lower critical bending moment, this increment 
is applied while the primary bending moment is allowed 
to increase above its lower critical value by an amount 
such that B is as small as possible; for the determi- 
nation of the upper critical bending moment, the incre- 
ment is applied while the primary bending moment is 
held constant at its upper critical value. 

The term C, called the initial torsional rigidity, is 
similarly defined as the ratio between an increment of 
torque and the corresponding increment of twist per 
unit length. For this purpose it is immaterial whether 
the primary bending moment remains constant or in- 
creases during the application of the torque increment, 
since in either case it will be shown, by an argument 
due to Neal,! that C remains constant at its elastic 
value irrespective of the amount of plastic flow which 
has occurred under the action of the primary bending 
moment. 

Whether we are concerned with the upper or lower 
critical bending moments, it is evident that the flexural 
rigidities A and B, defined as above, depend on the 
magnitude of the primary bending moment if the elastic 


limit has been exceeded, For this reason Eq. (2) has 
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LATERAL 


been arranged so that it defines the length of a beam of 
given Cross section and material which would buckle 
laterally at a given critical bending moment, rather 
than in the more usual form that defines the critical 


bending moment of a given beam. 


(3) THE OVERALL PRIMARY FLEXURAL RIGIDITY 

Before proceeding further, it is necessary to make 
certain assumptions about the stress-strain curve of 
the material (Fig. 1). Throughout this work it will be 
assumed, first, that the stress-strain curve is the same 
in tension and compression. Secondly, if the material 
is loaded to a stress beyond the elastic limit and then 
the stress is decreased slightly, the material will behave 
elastically with a modulus equal to the Young’s modu- 
lus E, as shown in Fig. 1. 

Two other moduli will be referred to in the subse- 
quent work—namely, the tangent modulus -7, which 
is the slope of the stress-strain curve at any point, and 
the secant modulus Fs, defined as the ratio of total 
stress to total strain. Obviously, these are both func- 
tions of the stress. 

We now suppose that the beam, of width 2) and depth 
2d, is subjected to a primary bending moment J, 
giving a primary curvature A,, and a stress and strain 
distribution as shown in Fig. 2. The maximum stress 
and strain are f and e, while their values at a distance 
y from the neutral axis are o and ¢, respectively. Then 
we have 

K, = e/d = e/y (3) 


The bending moment required to cause this distribu- 


tion of stress is given by 


d 
M, = 4b f oy dy 
0 


From Eq. (3) we have 
y = e(d/e), dy = de(d/e) (4) 


so that the above expression may be written 


4bd? e 
M, = —- oe de (5) 
e 0 


It is now convenient to introduce three nondimen- 
sional quantities, relating to the stress-strain curve, by 


9 ~e 
P=-— I o de 
Je Jo 


> ~ 


Q =i J oe de (6) 
4 e 
R= ok ae’ de 


It will be seen that the quantity P is the ratio of the 
area Qabc under the stress-strain curve to the area of 
the triangle Obc (Fig. 1). The quantity Q is the corre- 


the equations 
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sponding ratio of the moments of the areas about the 
a-axis, while R is the ratio of the moments of inertia of 
the areas about the c-axis. Thus, for a given stress- 
strain curve, P, Q, and R are functions of the maximum 
primary bending stress f. For a strain-hardening ma- 
terial, such as aluminium alloy, it is clear that, what- 
ever the value of /, the quantities P, Q, and R satisfy 
the inequalities 

l < x */s (7) 


Ig P<2 1480< 4/5 


If f is below the elastic limit of the material, then 
P=Q=R=1 (8) 


Using the notation (6), we may write Eq. (5) in the 
form 


M,-= (4/3)bd2fQ (9) 


On using Eq. (3), we can differentiate Eq. (5) with 
respect to K, to give 

one _ dM, oe 4 bd? f (3 — 20) 

dK, de 3 e 
It follows that if Q = */, the slope of the curve relating 
primary bending moment and curvature becomes zero 
and the beam becomes unstable under the action 
of M, alone. However, this can only happen if a suffi- 
ciently large drop in stress has occurred before the 
stress f was reached. It is possible, for example, in the 
case of annealed mild steel if the ratio of the upper and 
lower yield stresses exceeds 1.5. 

The overall primary flexural rigidity A is obtained 

from Eqs. (3) and (9); thus, 


A = M,/K, = (4/3)bd*(f/e)Q 
This may be written in the form 
A/A, = Q/S (10) 


where A, [= (4/3) Ebd*] is the elastic primary flexural 
rigidity and S is the ratio of the Young’s modulus to the 
secant modulus at the stress f; thus, 


S = Ee/f (11) 
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(4) THe INITIAL SECONDARY FLEXURAL 


We turn now to the determination of the initial 
secondary flexural rigidity B, defined as the ratio be- 
tween an increment of secondary bending moment 
dM, and the corresponding increment of secondary 
curvature dK». Whether the primary bending moment 
remains constant or increases during the application of 
dM», an increment of primary curvature dA, will be 
produced in addition to the secondary curvature dA». 


RIGIDITY 
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Thus, the neutral axis corresponding to the increment 
of total curvature will be inclined at an angle a to the 
horizontal (Fig. 3), where 


tan a = dK2/dk, (12) 


Two possibilities arise, depending on whether tan a > or 
< d/b, as shown in Figs. 3a and 3b, respectively. In 
either case it is clear that, in the shaded regions of the 
cross section, the material unloads—i.e., material that, 
under the action of , alone, was in tension receives a 
compressive strain increment or vice versa. In these 
regions, therefore, the ratio between the increments of 
stress and strain is equal to the Young’s modulus £, 
since it is assumed that the unloading process is elastic. 
Over the remainder of the section, however, the mate- 
rial is loaded further, and the ratio between the stress 
and strain increments is equal to the tangent modulus 
E, at the prevailing primary stress. 

Also, at any point in the cross section, the increment 
of strain is given by 


de = ydK, — xdKe (13) 


Consider first the case where tana 2 d/b. It will be 


seen from Eq. (12) and Fig. 3a that 


tan a = dK.2/dK, = d/\b (14) 


where \ < 1. Also, on using Eq. (13), the increment 


of secondary bending moment dJ7/, is given by 


d b *rby/d 
: dMz = f E f x(xdKy — ydK,)dx + Er / x(xdK. — ydK,) ax | dy 
- 0 Aby/d b 


After dividing throughout by dK», using Eq. (14), and performing the integration with respect to x, this expression 


becomes 


dM, ib d dy AB y3 
B = —= 2-—-3-—+— 
dK» 3 0 d d* 


The first of these two integrals may easily be evaluated. 


Eq. (3), the second integral can be integrated by parts. 


53 7 dy A3y3 
dy + Beit? +3 — —-— }<¢y 
; 3 Jo d d3 


Also, on substituting KE; = do/de, and y/d = e/e from 


Proceeding in this way, we finally obtain the expression 


B I 
= -1A33R + S — 4) — GA(P + S — 2) + &(S + 1) (15) 
B, 16S 
where P and R are given by Eqs. (6) and S by Eq. (11). The term B, is the elastic secondary flexural rigidity so 
that 
B, = (4/3) Eb*d (16) 
Consider now the case where tan a < d/b. From Eq. (12) and Fig. 3b we see that 
tan a = dK2/dk, = yd/b (17) 


where y < 1. 


The increment of secondary bending moment in this case is given by 


1 : b ydx/b : : : d 0 . , 
5 4M, =E / f x(xdK2 — ydK,)dy dx + f Er f x (xndK2 — ydK,)dx dy + 
ad 0 0 b 


b fd 
f f Eyx(xdKy — ydK,) dy dx 
0 ydx/b 
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LATERAL INSTABILITY 


This expression may be integrated by a procedure 
similar to that used in the derivation of Eq. (15) and 
finally yields the equation 
B l E 1 3/ 
= ) 
B, 16S J 


where fy is the stress corresponding to the strain ye and 


2 id 
Py => f o de 
Tr(ye) 0 


u (Sy + Ry — 2P») | (18) 


3 ye 
)y = : 1 9 
¢ rice ae de (19) 
4 ye 
Ry = = : 2d 
Fiye)® 4 oe’ de 
Sy = E(ye)/f, (20) 


It will be seen that the definitions of Py, Q,, Ry, and 
S, are exactly the same as those of P, Q, R, and S, ex- 


cept that f has been replaced by fy and e by ye. If fy 
is below the elastic limit of the material, we have 
Py=Q = Rh =S,= 1 (21) 


If fy exceeds the elastic limit, they are all greater than 
1. Further, if y = 1, we have P, = P, Qy = Q, Ry = 
R, Sy = S, and fy = f. Eq. (18) then gives the same 
value for B/B, as Eq. (15) with A = 1. Also, if f, and 
therefore fy, is below the elastic limit, both Eqs. (15) 
and (18) give B = B,. These considerations afford a 
check on the analysis. 

Having established Eqs. (15) and (18), which en- 
able the value of the initial secondary flexural rigidity 
to be determined for any inclination a of the neutral 
axis, we now consider the conditions at the lower and 
upper critical bending moments separately. 


(a) Conditions at the Lower Critical Bending Moment 


As discussed previously, the lower critical bending 
moment is the lowest value of the primary bending 
moment at which lateral deflections can occur with in- 
creasing bending moment. The inclination of the 
neutral axis is such that the initial secondary flexural 
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rigidity is as small as possible. Thus the problem is 
reduced to minimizing B, as given by Eqs. (15) and (18), 
with respect to \ or y as the case may be. 

Consider Eq. (18) first. Using Eqs. (21), we see that 
(Sy + Ry — 2Py7) = Oif fy is less than the elastic limit. 
Also, it is shown in Appendix (I) that, if fy exceeds the 
elastic limit, (Sy + Ry, — 2P,) 2 0. Hence, the 
smallest possible value of B from Eq. (18) is obtained 
when the unloading zones of Fig. 3b consist entirely of 
material in which no plastic flow has occurred and is 
given by 


(22) 


B/B, = 1/S 


It remains to show that Eq. (15) does not give a value 
less than this. On differentiating Eq. (15) with respect 
to A, we obtain the equation 


(d/dX) (B/B,) = —(3/16S) 
(1 — 2) BR+ S — 4) + 2QP+ S — 3R)] 


Since \ < 1, we have (1 — A”) 2 0. Also, since R 2 1 
and S 2 1, we have (3R + S — 4) 2 0. Itis shownin 
Appendix (II) that (2P + S — 3R) 2 0, so that, from 
the above equation, (d/d\) (B/B,) < 0. Hence, B de- 
creases as \ increases, and the smallest value of B ob- 
tainable from Eq. (15) is given by substituting \ = 1; 


thus, 


B/B, = (1/S) + (3/16S) (S + R — 2P) 


The analysis of Appendix (I) shows that (S + R — 2P) 
> 0, and it follows that the smallest value of B obtain- 
able from Eq. (15) is greater than that given by Eq. 


(33). 


(b) Conditions at the Upper Critical Bending Moment 


The upper critical bending moment is defined as the 
value at which, if no previous lateral deflection had 
occurred, the beam would buckle laterally under con- 
stant bending moment. The inclination a of the neu- 
tral axis is therefore given by the condition that, while 
the increment of secondary bending moment, dM, is 
applied, the primary bending moment is held constant. 


Iftan a 2 d/basin Fig. 3a, this condition is expressed by the equation 


d b Aby/d 
f E f y(ydK, — xdK2)dx + Ep f y(ydK, — xdKe a dy = 0 
0 Aby/d ; 


b 


Proceeding as in the derivation of Eq. (15), the above expression finally yields the following quadratic equation for 


Xd: 
3\°(33R + S — 4) — 8A(3B + S — 220) + 6(P + S — 2) =0 (23) 
If the material is entirely elastic, P = Q = R = S = 1 and Eq. (23) gives y = 0 as we would expect. Also, it 
is easy to show that Eq. (23) givesa root A < 1 only if 
(24) 


6P + 160+ 9R+S < 48 


If this inequality is satisfied the initial secondary flexural rigidity B is given by Eq. (15). 


An alternative expres- 


sion, obtained by eliminating the term in \? by means of Eq. (23), is 
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B/B, = (1/6S) [\(3 + S — 20) — 3(P +S —2)+3(S + 1)] (25) 


If tan a < d/bas in Fig. 3b, the condition that the primary bending moment remains unchanged is 


6 ydx/b d 0 
EB 7 f y(ydK, — xdKe)dy dx + i Er f y(ydK, — xdK»)dx dy + 
0 Jo ‘ Ua 


After integration, this equation becomes 
320 + (y7f1/f) (6Py — 16Q7 + 9Ry + Sy) = 48 (26) 


Eq. (26) defines the value of y, but the transcendenta’ 
nature of the equation would, in general, necessitate a 
solution by trial and error. Knowing the value of y, 
the initial secondary flexural rigidity is given by Eq. 
(18). 

If y = 1, Eq. (26) agrees exactly with Eq. (24). 
Also, if fy is below the elastic limit of the material, 
Eqs. (21) and (26) give simply Q = */:. This is im- 
possible for a strain-hardening material with a stress- 
strain curve of the type shown in Fig. 1, and, as shown 
previously, even if the shape of the stress-strain curve 
is such that Q can attain this value, the beam is then 
unstable under the action of MM, alone. Thus it follows 
that, if lateral buckling occurs under constant bending 
moment, the unloading zones shown shaded in Fig. 
3b must contain some material in which plastic flow has 
occurred under the action of the primary bending mo- 
ment. This conclusion is in agreement with Neal’s 
analysis for annealed mild steel beams. 


(5) THE INITIAL TORSIONAL RIGIDITY C 


It was shown by Neal,' both theoretically and experi- 
mentally, that, if a beam has yielded under the action 
of a uniform primary bending moment, the initial tor- 
sional rigidity C remains at its elastic value irrespective 
of the extent of the yield. Neal’s experiments were 
confined to mild steel bars, but his argument applies 
equally well to bars of strain-hardening material. For 
completeness this argument is given below. 

Hill, Lee and Tupper® showed that, if plastic flow 
has occurred in a material subjected to a system of 
stresses whose components are o;, oy, 02, Try, Tyz and 
Tzz, the increments of strain de,, de,, de,, dyz,, @Yy2z, and 
dy.z due to a small increment of stress obey the equa- 


tions 
2Gd(e, — 2) = d(o, — &) + (eo, — ade) 
2Gd(e, — é) = d(a, — &) + (0, — &)dg? (27) 
2Gd(e, — @) = d(a, — &) + (¢2 — s)deh 
Gdyzy = dty + trd¢} 
(28) 


Gdyyz = dtyz + Tyd¢e 
Gdy,2 = dtiz + alat 


where 


dé = (1/3) (de, + dey + de.) = da/3K (29) 


b d 
f f Ery(ydK, — xdK2)dy dx = 0) 
0 ydx/b 


and 
& = (1/3) (o, + oy + @;) 


In the above equations G is the modulus of rigidity and 
K is the bulk modulus of the material. 

In each of Eqs. (27) and (28), the first term on the 
right-hand side represents a recoverable elastic strain 
increment and the second term an irrecoverable plastic 
strain increment. The parameter d¢ defines the magni- 
tude of the plastic strain increment and can have a non- 
zero value only if the stress increment satisfies the con- 
dition for continuing plastic flow. 

Eqs. (27), (28), and (29) are based on three assump- 
tions that are approximately borne out in practice. 
The first is that the plastic strain increment corresponds 
to zero dilatation. Secondly, the principal axes of the 
plastic strain increment and of the total stress coincide; 
thirdly, the difference between any two of the three 
principal plastic strain increments is proportional to 
the difference between the corresponding total principal 
stresses. These assumptions [and therefore Eqs. (27), 
(28), and (29)] apply not only to a material exhibiting 
ideal plasticity, such as mild steel, but also to a strain- 
hardening material. 

Now, under the action of the primary bending mo- 
ment alone, it is obvious that the shear stresses 7,, and 
Tx are zero whatever the extent of the plastic flow 
that has taken place. It follows, therefore, from Eqs. 
(28) that, when the increment of torque is applied, the 
increments of shear strain dy,, and dy,, are purely 
The problem of determining the initial tor- 
i.e., the ratio between the increments 
is therefore gov- 


elastic. 
sional rigidity 
of torque and twist per unit length 
erned by exactly the same equations as the purely elastic 
problem, and it follows that the initial torsional rigidity 
C retains its elastic value. 

For a beam of narrow rectangular section, of width 
2b and depth 2d, Timoshenko’ gives the following ap- 
proximate expression for the elastic torsional rigidity: 


C = (16/3)Gb*d[1 — (0.63 b/d) ] (30) 


This equation gives accurate values for C if d/b > 2 ap- 
proximately. 

It will be seen that the above argument is valid 
whenever the shear stresses Tyz and 7,, due to the 
primary loading alone are zero over the part of the sec- 
tion in which plastic flow has occurred. This is true 
not only in the case of uniform primary bending mo- 
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ment but also when the primary loading takes the form 
of an eccentrically applied thrust. Neal also shows 
that it is nearly true for a beam of ideally plastic mate- 
rial when the primary bending moment varies linearly 
in the axial direction. In this case the constant shear L 
force is carried almost entirely by the elastic core of 


x i tne e +Or 
the beam. To prove this it is necessary to consider the 3 ‘ 
° . . ° fom = 
yield criterion, and the argument cannot easily be ex- — - 
J ¢ Cy a=lO 





tended to include strain-hardening materials. oak — 
nA=eo 


(6) MATERIAL REPRESENTATIVE OF ALUMINIUM ALLOY 











The stress-strain curve of an aluminium alloy can often 0-65 | 
be represented closely by an equation of the form i. ¥ Ee_ £[i+(¢)] 
os Gl Ne 

e = (¢/E) [1 + (¢/a0)"] (31) o-al 
where E is the Young’s modulus and gp is a constant 4 
having the dimensions of a stress. The index » is usu- 
ally fairly large (of the order of 10, say) so that, until “ST 
¢ approaches oo, the term (¢/ao)" is extremely small. s 
The effect of changing the value of m is shown in Fig. 
4, and it will be seen that the limit as » — © represents ° oe eee ne 2a sees ce 
an ideally plastic material with a yield stress equal to Et/s, 
vig FIG. 4 APPROXIMATE STRESS— STRAIN CURVES FOR 


Eq. (31) will now be used to obtain some idea of the ALUMINIUM ALLOYS 


difference between the upper and lower critical bending 
moments. It is evident from Eq. (11) that with a 


stress-strain curve of this form Also, from Eqs. (6) and (31), we find that 


a = l ( f 0)” (32) nN S _ ] 
+ (f/o P=1+-— | | 
from which we obtain the equation Sint 2 
f 5 — yi" 33) Q=1+ n R — 1] 4 (S — =| (35) 
= ol. — ? (ede = Ts) 
‘ ‘ S?Ln + 3 2n + 3 
Substitution of Eq. (33) into Eq. (31) then gives 7 (S — 1)2 (S — 1)3 
. R=1+-G]- += | 
e = (a/E) S(S — 1)" (34) ° o* E + 4 n+2 3n + 4 


Now, on using Eqs. (3), (10), and (30) and the fact that 4, = AK, we may write Eq. (2) in the form 


Md 2 (° B, d? i) [8 lity @& T" 7 
Lr \SB BP S 2 — (1.26b/d) b? (36) 


where, as before, A, = (4/3) Ebd’, B, = (4/3) Eb*d, and vis Poisson’s ratio. 
By virtue of Eq. (34), this becomes 


2Ed 2 , (0 B, d? ) TO l+y d? ] 
= = — S(S — 1)" — 1 — 1 37 
P ool T (¢ B bb? k 2 — (1.26b/d) b? (37) 


Also, from Eqs. (9) and (33), we have 
m = M,,/(4/3)bd?00 = Q(S — 1)'/" (38) 
It is easy to see from Eqs. (32) to (35) and the analysis of Section (4) that, for given values of b/d, v, and n, the 
right-hand sides of Eqs. (37) and (38) are functions of S only. Hence, the parameter m, which defines the critical 
bending moment, is a function of p, which defines the ratio between the depth and the length of the beam. 
In the elastic range, we have Q = S = 1, B = B,, ande = M,,/(4/3)Ebd*. Eq. (36) then gives the straight- 


line relationship, 
9 12 \ '/2 1 ]2 1/s 
, (5-2) | = ~~] (39) 
m a \b? y 2 — (1.26b/d) b? 


The value of p at the lower critical load is obtained by substituting B,/B = S from Eq. (22) into Eq. (37), 
while the value of p at the upper critical load is obtained by using Eqs. (18) or (25) as the case may be. 
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Numerical calculations have been made for a mate- 
rial in which vy = 0.3 and nm = 10. Two values of d/b 
namely, 6 and S—and the results 
are shown in Fig. 5. Using Eqs. (24) and (35), it was 
found that \ = 1 if S = 4.34. For any aluminium 
alloy such a value of S would represent a strain at 
The calcula- 


were considered 


the outer fibers well in excess of 0.015. 
tion was therefore confined to values of S < 4.34, in 
which case the neutral axis for secondary flexure at the 
upper critical bending moment is as shown in Fig. 3a. 
For each value of d/b in Fig. 5, the upper curve repre- 
sents the upper critical bending moment and the lower 
curve the lower critical bending moment. In the shaded 
regions betweer. the two curves, lateral deflections can 
occur with increasing load. The parameter m is a func- 
tion of S given by Eq. (38) and is independent of the 
value of d/b. The scale of Sis shown at the right-hand 
side of Fig. 5. 


It will be seen that no appreciable reduction in criti- 
cal bending moment below the elastic value occurs 
until S reaches the value 1.1 approximately. There- 
after, the reduction increases progressively as S in- 
creases. Also, it will be seen that, over a large part 
of the range considered, the difference between the 
upper and lower critical bending moments is of the 
order of 5 per cent. 
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REGIONS IN WHICH LATERAL DEFLECTIONS CAN OCCUR 





AERONAUTICAL 


SCIENCES—DECEMBER, 1952 
(7) COMPARISON WITH ANNEALED MILD STEEL anp 
DISCUSSION 


Neal's experimental work! on lateral buckling was 
confined to annealed mild steel bars, and it is interest- 
ing to compare the theoretical results just obtained for 
an aluminium alloy with those for annealed mild steel, 
The stress-strain curve in this case may be represented 
closely by Fig. 6, in which oo is the lower yield stress 
and yoo is the upper yield stress. If the primary bend- 
ing moment is of such a magnitude that the depth of 
the central elastic core is 2d, the strain at yield must 
be equal to pe as shown. It is then a simple matter to 
show that 


Py = 2y + p(u — 2), P=2+ p(u — 2) 


O=3rtp(u- 3), O=! 


{ { . | 
4) R= 3+ P (. ~ ;) 


On substituting these expressions, it can be shown that 
Eqs. (15), (18), (23), and (26) reduce to the correspond- 


w 


I| 


+? (. os ;) 


bo 


7. a 
Ry = 3 vy? + p’ (. — 
Sy = yu/ Pp, 3 =e p 


ing equations derived by Neal. 


If p is replaced by u/S, the above equations for P, Q, 
and R become 


P=2-+ (u/S) (u — 2) } 
Q = (3/2) + (w?/S?) [u — (3/2)] (40 
R = (4/3) + (u?/S*) [u — (4 3) 4 

Using Eqs. (40), numerical calculations, exactly 


similar to those previously described for aluminium 
alloy, were made for annealed mild steel in which up = 
1.4. Poisson’s ratio vy was again assumed to be equal 
to 0.3, and the same two values of d/b—namely, 6 and 
The results of these calculations 
lower 


S—were considered. 
are also shown in Fig. 5. In each case the 
curve represents the lower critical bending moment 
and the upper curve the upper critical bending moment. 
It should be noticed that, as soon as yield occurs at 
m = 1.4, the value of S jumps suddenly from 1.0 to 
1.4. Calculations were confined to values of S between 
1.4 and 2.4, since at the latter value the neutral axis 
for secondary flexure at the upper critical bending mo- 
ment is almost at the transition from Fig. 3a to Fig. 
3b. 

It will be seen that in this case the difference be- 
tween the upper and lower critical bending moments is 
much less than for the aluminium alloy, being of the 
order of 0.5 per cent. The value u = 1.4 was chose 
because this corresponds to the material used in Neal’s 
experiments. He found remarkably good agreement 
between the experimental critical bending moments 
and the theoretical upper critical bending moments. 
However, the difference of 0.5 per cent between the 
upper and lower critical bending moments corresponds 
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to approximately one increment of loading in Neal’s 
experiments, and it is therefore felt that the agreement 
that Neal obtained with the theoretical upper critical 
bending moment is not conclusive evidence that the 
jatter is closer to the collapse load than the theoretical 
lower critical bending moment. In fact, on the basis of 
the known behavior of struts under axial compression, 
it is felt that the reverse may often be true. 

It must be emphasized that this suggestion is intuitive 
and that no experimental evidence is available. How- 
ever, apart from this question, the use of the lower 
rather than the upper critical bending moment as a 
measure of the collapse load has two points in its favor. 
In the first place, it is far easier to calculate, since 
both the overall primary and the initial secondary 
flexural rigidities are given by simple expressions 
namely, A/A, = Q/S and B/B, = 1/S. Secondly, it 
is conservative, being somewhat lower than the upper 
course, 


c 


critical bending moment, although it is, of 
possible that small initial deviations from straightness 
which inevitably occur in practice may cause the 
beam to collapse prematurely before the lower critical 
bending moment is reached. In this event, however, 
the lower critical bending moment will certainly be in 
better agreement with the collapse load than the upper 
critical bending moment. 

In conclusion, one further point should be mentioned. 
The case of beams of rectangular section is not of great 
practical interest except insofar as it is a guide to the 
expected behavior of other sections. The problem of 
determining the lateral buckling loads for open sections, 
such as channels and I-beams, is complicated by the 
warping of the cross sections under torsion and the 
consequent increase in torsional rigidity when the tor- 
sion is nonuniform (see, for example, Timoshenko’). 
However, it is felt that a considerable simplification in 
the calculations may be effected by determining the 
lower critical load at which lateral deflections can 
start with increasing load rather than the upper critical 
load at which, if no previous lateral deflections had oc- 
curred, the beam would buckle under constant load. 


APPENDIX (I) 


It is required to show that (.S, + R, — 2P,) is never 
negative. Let 


wt 


= fy(ye)* (Sy + Ry — 2Py7) 


= E(ye)* + if ae? de — sre" ff o de 
0 0 


Now, if fy is less than the elastic limit of the material 


we have, using Eqs. (21), § = 0. Also, on differenti- 


ating the above expression, 


dé/d(ye) = 4(ye)n 


- 
n = E(ye)? — 2 f o de 
0 


where 


843 














pet----9 
Gi-*- 7" 7 ; 
' 1 | 
' i 1 
| ' | 
w i | l 
= i \ 
w | 
a ! | \ 
o ! ] 
pe Ve e STRAIN 


FIG. 6. STRESS- STRAIN CURVE FOR ANNEALEO MILO STEEL 





If f, is less than the elastic limit, this gives 7 = 0. Also, 


on differentiating, 
dn/d(ye) = 2(Eye — fy) = 2f+(Sy — 1) 


Since S, 2 1, we see that dn/d(ye) 2 0. Since n = 0 
when f/f, is less than the elastic limit (i.e., when Sy = 
1), it follows that » 2 0. Hence, d&/d(ye) 0, and, 
since = 0 when /, is less than the elastic limit, it fol- 
lows that £ 2 0. Hence, we conclude that S, + Ry — 


oF, 2 @. 
APPENDIX (IT) 


It is required to show that (2P + S — 3R) is never 


negative. Let 


¢ = fe(2P + S — 3R) 


= ie f ode + Ket — 2 ae’ de 
0 0 


If f is less than the elastic limit of the material, we 
have P = S = R= 1sothat¢=0. Also, on differen- 
tiating the above expression, 


d¢/de = 4fe?(P + S — 2) 


If f exceeds the elastic limit, we have P > 1 and S> 1 so 
that d¢/de > 0. It follows, since ¢ = 0 when / is less 
than the elastic limit, that ¢ 2 0. Hence, we conclude 


that 2P + S — 3R 2 0. 
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Creep Collapse of Viscoelastic Columns 
with Initial Curvatures 


HARRY H. HILTON? 
University of Illinois 


SUMMARY 


Collapse due to creep of a viscoelastic column is investigated 
by using the generalized viscoelastic stress-strain relationship. 
An expression is derived for the ultimate time of the viscoelastic 
column in terms of its arbitrary initial curvature, the applied 
load, the material properties, and the ultimate moment of the 
The results are then applied to two 
namely, the standard linear solid 


equivalent elastic column. 
specific viscoelastic bodies 
and the Maxwell body. 


SYMBOLS 


= = constants 

By = constants 

E = Young’s modulus, lbs. per sq.in. 
I = moment of inertia, in.4 

vi = length of column, in 


Ma = applied moment, in.lbs. 

M; = internal resisting moment, in.lbs 

Mur. = ultimate moment, in.lbs 

F = applied load, lbs. 

Pr = Euler load, lbs. 

als = load ratio P/Pr 

t = time, sec. 

w = lateral deflection, in. 

Wo = initial lateral deflection, in. 

x = coordinate, in. 

5o = maximum initial lateral deflection, in. 

€ = strain, in. per in. 

o = stress, lbs. per sq.in. 

Te = relaxation time of stress under conditions of con- 
stant strain, sec. 

To = relaxation time of strain under conditions of con- 


stant stress, sec. 


INTRODUCTION 


i pecanponed AND EXPERIMENTAL‘ EVIDENCE indi- 
cates that metals do not obey Hooke’s law at any 
level of stress, although in the lower stress ranges the 
strains are still completely recoverable. This anelastic 
behavior of the metal must be represented by relations 
between stresses, strains, and their time derivatives. 
It should be noted that the viscous behavior of metal 
grain boundaries become more and more prominent 
with increasing temperatures. 

The high speeds of modern airplanes have made it 
mandatory to consider the effects of elevated tempera- 
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tures (due to aerodynamic heating, power plants, ete. 
on the airplane structural components. Recently, 
several solutions®» * have been presented taking into 
account temperature effects by using viscoelastic stress- 
strain relations to describe the anelasticity of metals. 

The viscoelastic creep-buckling problem has been 
treated by Freudenthal,’ who obtains by successive ap- 
proximations a solution for the lateral deflection in terms 
of a power series involving the applied load and time 
and then defines that value of time which causes the 
series to diverge as the critical time. Kempner'® has 
found that the small deflection theory will not give a 
finite critical time for linear viscoelastic columns.  Al- 
though experiments on plastics’ have established that 
the buckling load is time-dependent, there does not 
seem to be an analytic solution available for estimat- 
ing the critical time. 

An attempt is made in this paper to overcome this 
difficulty by re-examination of the collapse criterion 
and finding instead the ultimate time, at the end of 
which the viscoelastic column will collapse due to creep. 


ANALYSIS 


Consider the viscoelastic column of Fig. 1, simply 
supported at both ends and with an initial curvature 
due to the initial lateral displacement wy) = wo (x). If 
the column is now compressed by a load P less than the 
ultimate load, the lateral deflection becomes w = 
w(x, tf). If the loading is so rapid that creep effects do 
not have time to materialize until the load P is reached, 
then the loading is essentially elastic and the lateral de- 
flection corresponding to P before the column begins to 
creep is given by the theory of elasticity’ as 


< on 


yoy | Sin (nwx/L) 
n=1 1 — (P n~) 


w(x, 0) = w(x) = 
(1 


where the 6, are the coefficients of the Fourier series 


wo(x) = >> bo, sin (nrx/L) (2 
n=1 
The linearized column theory is used, since, if the 
stresses are below the proportional limit as assumed, 
then small deflections result. 
In a uniaxial stress field and under conditions of ma- 
terial incompressibility, the general viscoelastic stress- 
strain relationship may be written in the form’: * 
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and E£ is Young’s modulus. The material constants a; 
and b; may be time-dependent but do not vary through- 
out the isotropic and homogeneous body unless a tem- 
perature gradient is present. 

In an elastic column of constant moment of inertia 
], the equation for the internal resisting moment is 


M, = EI [(0°w/dx*) — w"] (6) 


where the primes denote differentiation with respect to 
%. 

Using Alfrey’s elastic-viscoelastic analogy,' one may 
write the expression for the internal resisting moment 
of a viscoelastic column whose stress-strain relation is 


given by Eq. (3) as 
R {M,} = EIQ {(Q*w/dx*) — w"} (7) 
The applied moment at any section of the column is 
given by 
M, = —Pw (8) 
Up to the point of collapse, the column will be in a 
state of stable equilibrium such that the applied and 


internal moments will be equal. Substitution of Eq. 
(8) into Eq. (7) yields 











—k?R{w} = Q{(0°w/dx*) — wo" } (9) 
where 
k? = P/EI (10) 
A solution of Eq. (9) is in the form 

w(x, t) = g(t)-filx) + fo(x) (11) 
W.; + By substituting Eq. (11) into Eq. (9), the functions 
° f and g may be evaluated from the solution of the fol- 

lowing differential equations: 
Fic. 1. Column notation. {R 7 ,2Q} g(t) ~ (12) 
R {co} = EQ {e} (3) fil(x) + Ban2h(x) = 0 (13) 
where the differential operators are fo"(x) + (R2a0/bo)fo(x) = wo" (x) (14) 


i />yi i- i-1 
R = (0°'/dt') + (a-,0°"1/0t~") +... + ao (4) with the constants a, to be determined from the bound- 
Q = (b,07/dt) + (b;-,07-'/0P~") +... +50 (5) ary conditions. 


The deflection w(x, ¢) of Eq. (11) is, therefore, 


@ on . \ 4 
w(x, t) = p> $8 (au) [A, cos (ka,x) + B, sin (Ran»x)| + E —aeP* can | sin (nwx/L) ¢ (15) 

if bb * O and 
w(x, 1) = 3 {g(apt) [An cos (kanx) + B, sin (Rox)]} (16) 


n=1 
when b) = 0. The function g(a?) may be readily found from the solution of the differential Eq. (12), once a par- 


ticular stress-strain relationship is specified. 








One of the initial conditions is given by Eq. (1), and 
the others, if any, are dependent on the desired visco- 
elastic stress-strain relation. The boundary conditions 
for acolumn simply supported at both ends are 


w(0, ¢) = w(L, t) = 0 (17) 
These conditions determine the constants 

Ay ® (18) 

for bb) ¥ O 

B,, = [6on/g(0)] $1/[1 — (P*/n?)] - 

1/[1 — (aoP*/bon?)|} (19) 

for bo = 0 
By, = [6on/g(0)] {1/[1 — (P*/n?)]} (20) 
a, = nw/kL (21) 


If time derivatives of higher than the first order are 
present in the chosen stress-strain law, then the addi 
tional constants in the function g(a@,f) are determined 
from the additional initial conditions associated with 
the chosen viscoelastic body. 

A column will collapse as soon as a state of unstable 
equilibrium is reached—i.e., whenever the configur- 
ation of the deflected column reaches the shape where 
the internal resisting moment is less than the applied 
moment. 
elastic column, the type of loading present may be 


Once creep begins in the deflected visco- 
considered as combined bending and compression. 
In the elastic case, for P < Py, the ultimate moment 
for such a loading may, be estimated from Shanley’s 
interaction curve,* such that 


Mune 


(Ia,/c) [1 — (P/Ac-,)* |" (22) 


P 2. 


) 


with B, given by Eq. (19). 


If the initial curvature is in the form of a single sine 
wave, then Eq. (25) simplifies to 


7.) (1/P*r,)] 


fat.[ C1 € 
[((1/P*) — 1] In (Myy,./PB:) (26) 


)- 


For the Maxwell viscoelastic body, Eqs. (25) and 


(26) are equally applicable, except that r, is infinite and 
B, is determined from Eq. (20). 
coelastic column all the quantities in the last two equa 


For any given vis- 
tions are known, and the ultimate time f,, corre 
sponding to a load P or vice versa may be readily 


found. 


Similar expressions can, of course, be derived from 
Eqs. (15) and (16) for different boundary conditions. 
It must be remembered that the ultimate moment, the 
relaxation times, and Young’s modulus are functions 
of the temperature and that, consequently, the ulti 
inate time is also temperature-dependent. 


S46 


(B, exp ftur[U/7.) — (a 


where o, and o,, are, respectively, the bending modulys 
of rupture and the compressive yield stress of the 
column of cross-sectional area A and with ¢ represent. 
ing the distance to the point furthest removed from the 
neutral axis. The exponents p and gq are character. 
istics of the particular material. 

Applying this collapse principle to the viscoelastic 
column, one sees that, at some time /,,, the deflection 
will be of such magnitude that the maximum applied 
moment (at x = L/2, for the simply supported column 
will reach the ultimate value described above. 


r 2, (g(a, alt) PB, 


n=1 


M an (23 


then permits the determination of the ultimate time 
t 
to creep after the application of the load P. 


1.e., time necessary for the column to collapse dui 


ult 


DISCUSSION 


For illustrative purposes, the results will be applied 
to two relatively simple viscoelastic media, which 
describe adequately the behavior of some metals and 
plastics. They are (a) the standard solid® and (b) the 
Maxwell body.’ 

In the case of the standard linear solid, the coeffi 
cients of the differential operators QO and R are 


a = 1/7. @ 0 (7 > 1)) 
by l/r,, 0d, | (24 
b, oS > ) 
where 7, and 7, are the relaxation times under condi 
tions of constant strain and constant stress, respec- 
tively. Eq. (23) now reduces to 
2/P*r,]/[(n?/P*) — 1]}) M un (25 
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A Note on the Contribution of Configuration 
Asymmetries to the Free-Flight Motion of 
Missiles 


John D. Nicolaides and Leonard C. MacAlister 
Aerodynamicists, Free Flight Aerodynamics Branch, Ballistic 
Research Laboratories, Aberdeen Proving Ground, Md. 


August 15, 1952 


ae FREE-FLIGHT MOTION of missiles having slight configur 
ational asymmetries* (i.e., control surface deflection, bent or 
damaged wings and/or fins, manufacturing inaccuracies, bent or 
) has been the subject of both theoretical 
The results of 


damaged bodies, ete 
and experimental investigation! at this activity. 
the theoretical investigation, employing a coordinate system fixed 
in space, indicate that, for the case of a missile flying at constant 
axial velocity and constant rolling velocity, the free-flight pitching 
and yawing motion is 7vricyclic [i.e., the motion may be con 
sidered to be comprised of three rotating vectors that are damp- 
ing or expanding (see Fig. 1)] and that the free-flight transverse 
displacement of the mass center is also Tricyclic plus a linear term 
and a constant. When the asymmetries are removed, the mo 
tions reduce to the familiar Epicyclic case.?~7 

In addition, the theoretical results indicate that one of the 
effects of configurational asymmetry is to introduce a resonance 
phenomenon when the rolling frequency of the missile approaches 
coincidence with one of the natural frequencies of pitch and yaw 
nutation rate This resonance phenomenon may also be ob 
tained from the results of reference 8, which employs a moving 
coordinate system, by noting that 


> 


Cu — 8P,) (Xe — tPe (1) 


(ac — bd) = Ayre 
where \; and d» are given by Eq. (19), reference 8, corrected for 
typographical errors as 


nN = [-—(a + c — eb)/2| + 


9 


(1/2) Via +e-— eb)? — 4(ac — bd (2) 


and where \; and \» are equivalent values of \; and 2 in a coordi 


nate system not rolling with the missile. 


Eq. (20), reference 8, may thus be written as 


(cs = bt) 
3: Ce , A A ‘ 1 
Rr + lr — P:)) [RA2 + 1d2 — P;)) 
from which it is seen that if RA; = 0 and P, + J); then n, > 
However, inspection of Eq. (1), Eq. (2), and the quantities as 
given in reference 8 indicates that it is improbable that RA; = 0. 


Nevertheless, the value may be low enough to yield large values of 
n: Near resonance; using the steady-state trim, zp,—g 25a meas 

* The contributions of the asymmetries to the aerodynamic force and 
moment system are assumed to be linear, and, when they are removed, it is 
required that the basic configuration have trigonal or greater rotational sym 


metry and mirror symmetry 
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Tricyclic pitching and yawing motion 


ure, magnifications of 10 — 20 by the resonance phenomenon 
appear usual and values of 50 —~ 100 are not unreasonable. Need 
less to say, the recognition and calculation of the effects of this 
phenomenon on stability play an important role in the contem 
porary design and testing of winged and /or finned missiles 

A preliminary program was undertaken in 1949 in the Aber 
deen Spark Photographic Range to determine the ability of the 
Tricyclic theory to represent accurate experimental data The 
results! indicated that the pitching and yawing motions could be 
represented by the theory to within the accuracy of measurement 


in the Range (linear measurement, 10~? ft angular measure 


ment, 7 min.; and time measurement, 10~® sec 
An experimental program for a detailed investigation of the 


resonance region is now in progress. 
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Prediction of Yawing Stability Characteristics 
of Airplanes During Catapulting—Errata and 
Additional Remarks 


Henry J. Kelley 

Research Department, Grumman Aircraft Engineering Corporation, 
Bethpage, N. Y. 

August 20, 1952 


— WRITER WISHES TO STATE correctly Eq. (4) of his article 
of the above title in the August, 1952, issue of the JouRNAL 


and its linearized counterpart Eq.(8). These should read: 


—(Iy + Mc?)v — McX cos W — F(b + c) sin (@+ Vv) + 


Mc¥ sin ¥ — P,\l =0 (4) 
—(ly + Mc)e — Mcé — Fb +c) (ye + ¥) + 
MYoy —Pl=0 (8) 


The additional terms appearing in the original do not affect the 
subsequent results in any way, since they were not carried through 
to the working equations. 

Of considerable interest is a method for estimating tire stiffness 
devised by K. R. Thorson, of Boeing,! which has recently come to 
the writer’s attention. This may be regarded as a decisive step 
in making possible reliable prediction of catapulting stability 


characteristics. 


REFERENCE 


1 Thorson, K. R., A Rational Method For Predicting Tire Cornering Forces 
and Lateral Stiffness, Boeing Report D-11719, March 28, 1951. 


On the Noncirculatory Flow About a Two- 
Dimensional Airfoil at Subsonic Speeds 


Bernard Mazelsky 

ae <> area Division, Lockheed Aircraft Corporation, Burbank, 
alif. 

August 11, 1952 
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tack, the forces and moments due to apparent or virtual magg 
flow are impulsive in nature and can be defined conveniently by 
the Dirac delta function. 

For the case of compressible flow, the lift and moment due to g 
sudden change in angle of attack are no longer impulsive but are 
finite and time-dependent. The association of the concept of 
apparent or virtual mass to the noncirculatory flow for the com. 
pressible case is questionable and requires examination. An jp. 
vestigation of the applicability of the apparent mass concept for 
subsonic flow was recently made by Miles.! The calculated re. 
sults, based on the acoustic approximation, were made for a circy- 
lar cylinder subjected to a suddenly applied force and a suddenly 
applied velocity that is initially at rest. The results are ap. 
plicable for extremely small Mach Numbers. 
in his paper that the concept of virtual mass applied only to the 
steady state and not the transient when the flow is compressible 

In the present note, the limitations of the virtual mass concept 


Miles indicated 


noted by Miles are illustrated for the case of a two-dimensional 
airfoil (flat plate) in a subsonic compressible flow that is initially 
flying at a relatively high subsonic Mach Number (0.5 < M < 
0.7). The noncirculatory lift and moment on a harmonically os- 
cillating airfoil in a subsonic compressible flow can be obtained 
in closed form in terms of Mathieu functions. Details of this 
evaluation for the cases of sinking and pitching motion are given in 
several papers (for example, see reference 2). Complete tables 
of the Mathieu functions required for this calculation have been 
published recently. Results for the lift due to sinking velocity 
are given in Table 1 as a function of the reduced frequency, k, for 
three Mach Numbers. The results are given in terms of a com- 
plex lift derivative, F,(k) + iG.“(k), which is defined by the 
following equation for lift: 


L® = rpcV? e*ikh [F.O(R) + iG.(R)] (1 


where k is the reduced frequency, wc/2V, h is the amplitude of 
vertical displacement in half-chords, and s is the distance traveled 
in half-chords. The superscript (1) denotes the noncirculatory 
portion of the flow. 

With the aid of the following reciprocal relations, the noncircu- 
latory lift derivative, ki)(s), due to a sudden acquisition of verti- 
cal velocity can be calculated from the results for the harmonically 
oscillating case given in Table 1: 


2 © F(R) sin ks 
7 CALCULATION OF THE TRANSIENT AERODYNAMIC PROPER- kis) = dk (s>0O) 
TIES of an airfoil requires consideration of a circulatory and a FF k 
noncirculatory flow. The present note is concerned with only Po (2 
the noncirculatory flow at subsonic speeds. For the incompres- : = 
. ° . € * ¢a) ~ » . 
sible case, this flow is commonly referred to as the apparent or iis) ws = G.(k) cos ks oe 
virtual mass of the airfoil. For a sudden change in angle of at- ™ Jo k 
TABLE 1 
Complex Noncirculatory Lift Derivatives Due to Oscillatory Sinking Velocity 
- -M = 0.5- ~— —M = 0.6- ~~ — M = 0.7 - 
k F(R) G(R) k F.(k) G(R) k F(R) G(R) 
0 0 0 0 0 0 0 0 O 
0.67083 0.0378 0.43824 0.47703 0.02894 0.3564 0.32582 0.0220 0.2532 
1.06059 0.1801 0.7371 0.76603 0.1389 0.5719 0.51518 0.1030 0.4270 
1.50000 0. 5662 0.9907 1.06667 0.4244 0.7521 0.72857 0.3146 0.5661 
1.83711 0.9218 0.9512 1.30640 0.6822 0.7220 0.89231 0.4981 0.5433 
2.12133 1.1096 0.7865 1.50850 0.8116 0.6045 1.03036 0.5843 0.4607 
2.59809 1. 1646 0.5427 1.84752 0.8305 0.4418 1.26193 0.5847 0.3596 
3.35409 1.0817 0.4985 2.38514 0.7815 0.4713 1.62914 0.5536 0.4267 
3.67422 1.1215 0.5474 2.61278 0.8330 0.5297 1.78463 0.6063 0.4845 
4.24263 1.3000 0.5471 3.01700 1.0148 0.5343 2.06071 0.7680 0.4908 
4.74339 1.3782 0.4104 3.37310 1.1158 0.4083 2.30394 0.8452 0.3985 
5.61249 1.3578 0.1655 3.99109 1.0904 0.2825 2.72606 0.8207 0.3330 
6. 70802 1.3279 0.1712 4.77028 1.1302 0.2655 3.25828 0.9287 0.3478 
8.21588 1.2821 0.1025 5.84238 1.1294 0.1168 3.99055 0.9641 0.2225 
9. 48684 1.2338 0.08424 6.74619 1.0863 0.0919 4.60789 0.9861 0.1758 
11.61894 1.2565 0.1274 8 . 26236 1.0376 0.0638 5.64349 0.9493 0.0936 
13.41642 1.2714 0.1361 9.54055 1.0322 0.1371 6.51654 0.9223 0.0537 
15. 00000 1.3014 0.08362 10. 66667 1.9326 0.1008 7.28571 0.8966 0.0594 
© 4/7 0 x 1/0.37 O co 2/0.79 0 
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FIG.I. NON-CIRCULATORY LIFT COEFFICIENT DUE 
TO A SUDDEN CHANGE IN SINKING VELOCITY 


where k,‘(s) is defined such that 


L® = arpcV%(dh/ds)ky(s) (3) 


The results for &,(s), which were calculated by evaluating Eq. 

(2) numerically, are plotted in Fig. 1 for three subsonic Mach 

Numbers. The values of &,“(s) at s = O [which is lim F.()(k)] 
k—> = 


were computed by the following equation: 


k,(0) = 2/7M (4) 


Also shown on Fig. 1 is a graphical representation of k"(s) func- 
tion at s = Owhen M = 0. 

It is of interest to note that the general shapes of the curves in 
Fig. 1 are similar to the result obtained by Miles for a circular 


Note on the Determination of Bending 
Frequencies by Rayleigh’s Method 


F. Niedenfuhr 
North American Aviation, Inc., Downey, Calif. 
August 25, 1952 


M**" STRUCTURES THAT ARE ANALYZED by simple beam theory 
have mass and stiffness distributions that are not constant. 
In particular, one frequently finds systems like the cantilever 
beam shown in Fig. 1 which is heavier and stiffer near the reot 
than near the tip. 

A simple Rayleigh analysis to find the first natural frequency 
of the beam shown in Fig. 1 might proceed as follows: 

First, a mode shape is selected which satisfies the boundary 
for instance, the fourth-degree ex- 


conditions of the problem 


cylinder kubjected to a sudden velocity (see Fig. 5, reference 1). 
Although direct comparisons cannot be made because of the 
differences in body shape and initial conditions, qualitatively, the 
conclusions indicated by Miles also appear applicable to the re- 
sults obtained in the present investigation. 


REFERENCES 


1 Miles, John W., On Virtual Mass and Transient Motion in Subsonic Com- 
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2 Reissner, Eric, On the Application of Mathieu Functions in the Theory of 
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pression 


x\? 1 /x Ll sx} 
y= ( ) 2—- T y] (1) 
‘ \/ 3\l 3 
is seen to satisfy the boundary conditions 
y (0) = y(0) = 0] 
- : (2 
y"(l) = y'"() = of ) 


Next, the beam is considered to be divided into a number of sec- 
tions, and the mass, stiffness, deflection, and curvature (y’’) are 
calculated at each section. These quantities are shown in Table 
1, Columns 3, 4, 5, and 6. 
equal to the lateral area of that section, and the stiffness is taken 


Here, the mass of a section is taken 


equal to the cube of the depth. 
The calculations indicated in Columns 7 and 8 of Table 1 are 
then completed, and the frequency is found by means of Ray- 


leigh’s formula 
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w* = ZEI(y’’)?/Zm(y)? (3) 


The result is shown at the foot of Table 1. 

The validity of this process and the accuracy of the result de- 
pend critically on the accuracy with which one is able to pick the 
mode shape, y. As is well known, however, small errors in a 
curve may result in large errors in the higher derivatives of that 
curve, so that the frequency as given by Eq. (3) is in error by an 
amount proportional to the error of the second derivative of the 
mode shape that waschosen. This error is generally much greater 
than the error in the mode shape itself. The state of affairs can 
be seen by inspection of Table 1, where large values of stiffness 
get multiplied by large values of curvature, and the numerator of 
Eq. (3) is almost completely determined by what happens at the 
root end of the beam. 

An alternate to the above procedure is to construct a deflection 
curve starting from an assumed curvature. A simple method of 
guessing a reasonable second derivative to the deflection curve is 
as follows: 

Begin with the Bernoulli-Euler Relation, which states 









































































































































curvature = y’ = M/EI (4) 
TABLE 1 ‘| 
‘ 2 3 a Ss 6 7 8 
ni] x m el y" y Ewyy® | my)? 
10-4 
‘ 5S 19 6859 0361 004835 |89.387174 | 000 044 
2 Ls 17 4913 0289 -.040 669 |41.033867 | .002 812 
3/25 LS 3315 0225 105 469 |17685937 | .c166686 
4 35 3 2.197 .0169 192 835 | 6274852 | .048341 
5 | 45 uw 1.331 0121 297169 | 1.948 717 | .097140 
6 SS 09 0.129 .00568 326035 | .235193 | .09oSs68) 
7 65 0.7 0,343 .0049 538335 | 082354 | .202685 
8 | 7s 0.5 0.125 0025 .667969 | -0°7813 | 223001 
9 | 85s 04 0.064 .0009 .B001e9 | .000518 | .256108 
10 |} 95 | 0.4 0.064 .0001 .933336 | .000006 | .398.446 
u* = 150.036.35' 465“ 5 w.ito y I tn 1.291212 
Tassie 2 
' 2 3 4 5 e 7 8 
n x m EI y” y Er(y") | miy)* 
‘ 5 1.9 6.859 13748 -018 129659 | o00 616 
2 Ss 1.7 4913 16853 178 139 516 053 862 
3/25 Ss 3.315 21037 491 449 383 | -26162) 
4 35 13 2.197 26991 1.014 160049 |!.336655 
Ss 4s i 1.33 35687 1Bi6 169 384 | 3.627 642 
6 55 09 0.729 49383 2.984 177780 | 8.013830 
7 165 0.7 0.343 -699]0 4.653 167 923 [15.155 286 
8 75 os 0125 1.12000 7043 \S6800 |24.801925 
> 85 04 0.064 \.000 00 10.352 064000 |42.865560 
10 | 9S 0.4 0.064 23410 | 14.454 003525 |93,567246 
at: Re Or } W=.0857 s 2 (MBO /1794.424242 
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where ./ is the bending moment. Now, E/ is known everywhere, 
and J is of nearly the same shape no matter what the shape of 
the deflection curve is, for a particular set of boundary condi 
tions. Thus, a simple shape should be picked for / and y"’ calcu- 
lated from Eq. (4). 


and the Rayleigh method can proceed as usual. 


A double integration then determines y, 


Fig. 2 shows an assumed bending moment and the resulting 
curvature for the beam of Fig. 1. This curvature was numeri 
cally integrated twice to obtain the values of y shown in Table 2 
The frequency is 


y”’ in Table 2 was taken directly from Fig. 2. 
shown at the foot of Table 2. Note that the boundary condi 
tions at x = 10 have been taken care of by choosing WM = dM/dx 
= Oat x = 10 and that the boundary conditions at x = 0 have 
been taken care of by making the constants of integration zero. 

The frequency obtained by the first method is 28 per cent higher 
than that obtained by the second method. That the lower value 
is more correct can be argued as follows: 

The calculations are the same except for the choice of deflection 
function, and Rayleigh’s theorem states in effect that, other things 
being equal, the mode shape that yields the least frequency is 
closest to correct and that, whatever mode shape is chosen, the 
frequency calculated from it will be too large 

The method of picking the bending moment as proposed here 
compares fairly well with the first method with regard to time 
required for computation. Approximately 20 per cent more 
time was spent preparing Table 2 than Table 1, although the 
time required for Table 1 did not include the determination of the 


numbers for mass and stiffness 
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Simplified Laminar Boundary-Layer 
Calculations for Bodies of Revolution and for 
Yawed Wings—Errata 


N. Rott and L. F. Crabtree 
Cornell University, Ithaca, N.Y. 
September 15, 1951 


UR ATTENTION WAS CALLED to the following errors in our 
O “Simplified Laminar Boundary-Layer Calculations for 
Bodies of Revolution and Yawed Wings” (Journal of the Aero 
nautical Sciences, Vol. 19, No. 8, pp. 553-565, August, 1952). 


” 


In Eq. (5.3), the exponent of v should read 0.2, instead of —0.2, 
and the exponent of U’ (outside the bracket) must be changed 


from —3 to —3.4. The same corrections apply to Eq. (5.4) 


Note on the Two-Dimensional Losses in 
Cascades 


A. Charwat 
Propulsion Research Corporation, Inglewood, Caiif. 
August 25, 1952 


ITH REFERENCE TO Eg. (10) of the paper on ‘‘Two-Dimen 
pe Losses in Turbine Blades’’ by MacGregor,! relating 
the energy loss factor through a cascade to boundary-layer 
parameters, it is advantageous to redefine the ‘‘energy thickness” 


as follows: 


1 1 6 
me puy* gi = a* | (1493 u3)dy 
~ - J0 


This form follows directly from the consideration of the deficiency 
of kinetic energy of a stream of profile u referred to one of con 
stant velocity mp. In particular, if up is the mean outflow velocity 
and the integration is carried over the pitch of the cascade, the 
continuity equation 


~, - 
{ udy = f udy 
0 0 


« . 


allows it to be rewritten as 
as si 
a u \2 
= l dy (1) 


#1 
Uo 


0 1g 
where ¢; represents the energy lost per unit breadth and unit 


time 
Under these conditions the blade loss factor can be written 


directly as 
, energy lost/unit volume 4 ( l ) 
" energy available/unit volume Z energy available , 
(ae lost /unit time ) _ MI 5) 
volume flow /unit time r Ugl ue 
Eq. 2 above and 


using the notation of the mentioned article. 
Eq. (10) of reference 1 are shown to be identical expanding ¢): 


tu n \? u u 
1 = 1 dy = oa = dy + 
J 0 to uo Q Mo uo 
‘fu \? u 
1--— dy =@0+~¢ 
o \to Uo 
hence, 
9/6 =1+(1/K) =(K +1)/K (3) 


and using the definition ug = to cos ¥. 
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If the energy lost is taken to equal the work done by the shear 
forces, then 
I 


= ply ¢; = Up rT) dx (4 ) 


J 0 


Imagining that an ideal flow up to the trailing edge is changed 
at constant pressure—that is, constant free-stream velocity—to 
a flow with a momentum corresponding to that of the actual flow, 
we can use von Karman’s momentum equation to compute the 
required work expenditure to perform the change between this 
(fictitious) condition (subscript 1) and the actual condition 
(subscript 2). 


*) 


9 
i 

uo Tm dx = u ply? = pity 
l ss 1 dx 


since 0, = 0. Hence, from Eq. (4) 


¢1/0 2 
The above corresponds to the value of A 1, and the loss factor 
w becomes 
WwW = 20uo/tial (5) 

This relation has been implied before [reference 2, Eq. (2.12)] on 
the basis of the same assumptions, of course, but by use of a dif- 
ferent argument. 

The fictitious flow considered, obviously neglects the dissipa- 
tion of energy due to gradients in the boundary layer, as implied 
by the value of unity for K 


' t 
u(u u) dy 
uy? J 9 


K =1 7 
l 9 
: u?(uo — u)dy 
No? 0 
which, expanded, leads to 
*! 
u(u uy)? dy = 0 
0 
and since the integrand cannot be negative, uo wu. Tomakean 


order-of-magnitude estimation of the dissipation in the two 
boundary layers of the blade passage, we assume the values for a 
turbulent boundary layer on a flat plate to represent the average 
conditions between the thin layer on the pressure side and the 
thick one on the suction side of the blade passage. Thus, the 
total dissipation per unit time is 


ae 
> 


*e ° ou \? 
y == 2 ob dy dx: o= “ > 
0 0 oy 


Using the '/; power law* 


Uo? 
y = 2 / dx a 
J 0 ‘ 


bp = —0.77uouR., 


Using farther 


ix) = 


Now, if the deficiency in energy of the flow equals the work done 
by surface friction plus the dissipation into heat then, 


(1 2) pu 41 = ply’O 4 y 
1 K+ 1 us i ; 
= : = 2 1.54 - R, 
0 K 0 
For average conditions of R,. = 6.5 X 10° and a mean value of 
C = 3in., 6 = 0.0125 in., 
¢gi/@ = 1.9918; K = 1.009 


The results are compared with predictions for the various 


values of K on Fig. 1. The last computation yields, of course, a 
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Fic. 1. Comparison of calculated and measured loss factor (+) 
K = 1.280, reference 1 (@) K = 1.00, and Eq. (6). 
value of K which is not constant. However, that correction is 
small, indicating that the dissipation in the boundary layer and, 
consequently, that in the blade wake is small and that the much 
simpler relation (5) can be used to a fairly high accuracy in esti- 
mating two-dimensional cascade losses. 
REFERENCES 
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On the Stability of Laminar Flame Fronts* 
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¥ REFERENCES 1, 2, AND 3, the stability of a laminar flame front 
to two-dimensional disturbances that are time-variable is 
investigated. In references 1 and 3 the flame front is considered 
as a discontinuity in velocity, whereas in reference 2 a continuous 
distribution of velocity in the direction of flow is considered. 
In all three works, the effect of energy dissipation is neglected. 
A somewhat more general consideration of the problem is herein 
presented. 

The equations representing the flow of a viscous fluid may be 
given as follows: 
Continuity: (Op/dt) + (pux),~. = O 


p(du,;/dt) = —pye + (us. be), e+ (mee ie — 
2/3) (wuz, ki 


Force: 


9 
Energy: p(dE/dt) = (KT.x).,. — (> + 3 Hi .) rT 


M( UK, Me, t+ UK, it, &) 


* This work was carried out under O.N.R. Contract NONR-656(01). 
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where the variables are defined as follows: 


x; = rectangular cartesian space coordinates 
t = time coordinate 

u; = velocity component in 7 direction 

p = density of fluid 

p = pressure 

T = temperature 

E = internal energy 

m = coefficient of viscosity (absolute) 

K = conductivity (thermal) 


If it is assumed that the flow consists of a parallel steady state 
portion plus a small, linear, three-dimensional disturbance and jf 
the disturbance is analyzed into its Fourier components, the flow 
variables and a single disturbance component can then be repre- 
sented as 

u, = U(x) + u(x) exp[i(axe + Bx; — vt)] 
Us = W(X) exp[i(axe. + Bx; — vt)] 
us = w(x) exp[t#(ax2 + Bx; — vt)] 


= p(x) + p(x) exp[i(axe + Bx; — vt)] 


p 
b = p(x) + p(x) exp[i(ax, + Bx; — rt)] 
T = T(x) + T(x) exp[i(ax, + Bx; — vt)] 
E = E(x) + E(x) exp[i(ax, + Bx; — vt)] 
mw = p(x) + u(x:) exp[i(ax. + Bx; — vt)] 


= K(x) + K(x) exp [i(ax2. + Bx; — vt)] 


The equations for the disturbance component may be obtained 
by substituting the above quantities into the equations of motion, 
subtracting the steady-state terms as satisfying themselves, and 
discarding the terms nonlinear in the disturbance. 

It should be noted that v and @ appear symmetrically with 
wand B. 

Continuity: —ivp + plu’ + iav + Bw] + 
pu + pU’ + p'U =0 


Force: p{[(U' — iv)u + Uu'|] + pUU' 


4 
—p'’ +h E u” — (a? + B?)u + : av’ + ow’) | + 


4 , =; 4 acd tye 
a’ F u ~s i(av + Bw) + 3 (ul + p’U') 


p(—ivy + Uv’) = —iap + cE — (a? + B*)v + 


t ae 

: (u’ + iav + ise) + p’(v'’ + tau) — 3 iapl 
p(—ivw + Uw’) = —iBp + Aw" — (a? + B2)w + 

ip , ; } . ae 

3 (u’ + iav + iBw) | + p(w’ + Bu) — : iBul 


e 


Energy: pUE’ + p|—ivE + UE' + uk’) = RT’ + 
RIT” — (a? + BT] + K'T’ + KT" - 


} 4 o 
(> he 3 4v’) (u’ + iav + iBw) + 4fU'u’ + 3 pu( U')*? — pl 


If the second force equation is multiplied by a and the third force 
equation by 8 and if both equations are then summed, the en- 
semble will then be in two-dimensional form if the following group- 
ing is performed similar to that of reference 4: 
av + Bw = ai 
a? + B? = a! 
The conclusions from the above are that, if a minimum Rey- 
nolds Number for stability of laminar flame fronts (based on 
flame front thickness) exists, it is the same for two- and three- 


(Concluded on page 864) 
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(Concluded from page 852) 


dimensional disturbances. Hence, it is sufficient to consider 
only two-dimensional disturbances. 

If the steady-state velocities involved in the problem are low 
compared with the velocity of sound, it is considered that the 
disturbance is isochoric and, therefore (for the two-dimensional 
disturbance), 

u’ + tav = O 
The two-dimensional disturbance equations therefore become 


Continuity: —ivp + p'u + (pU)’ = 0 
u’ + tav = O 


Force: p((U' — tv)u + Uu'| + pUL’ = 
—p' + plu” — atu) + p’[2u’] + 


f yyre 
= [wU" + wl’ 


p|—iw + Uv’) = —iap + plv” — av] + 


9 


a’ lv’ + iau) — = tap’ 
o 


For the isochoric disturbance, the energy equation need not be 
considered. However, it is necessary to define an equation of 
state for uw. If it is assumed that the disturbance does not affect 
the flame structure, then it follows that 7 = f(p) and wp = p(da 
dp). The boundary conditions for the disturbance are that 


u(+o) = Vor that 


u'(—«o)/ul—-~) =a; u'(+o)/u(4+o) = 


= 


Since two boundary conditions fully determine the problem, 
the general solution consists of a combination of two linearly 
independent solutions. 

u= Cim + Cotte 


The formal statement of the eigenvalue problem in dimension- 


less form can therefore be given as follows: 


u,'(— 2) —au(—@) 

, : 
‘ U2 (—-2) — ale —o) 
F(a, v, R) = 


uy’ ( +o) + an,(+ 


ux'(+0) + ato(+~) 


When detailed laminar, steady-state flame data become avail- 


able, calculations can be performed. 
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